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1 . 


INTRODUCTION 


This report describes the results of a three phase 
program of atmospheric studies related to aerospace activ- 
ities and remote sensing technology performed by Atmospher- 
ic and Environmental Research, Inc. (AER) under the spon- 
sorship of NASA/Langley Research Center (contract no. NAS1- 
15943) . The period of performance for this work is August 
1979 - May 1980. Parallel efforts were undertaken to 
investigate: (a) the sensitivity of one-dimensional photo- 

chemical model simulations of projected supersonic aircraft 
operations to chemical rate constant data and parameteriza- 
tion of vertical eddy diffusion, (b) the feasibility of the 
development of a two-dimensional modeling capability based 
on the Generalized Lagrangian Mean (GLM) methodology, and 
(c) the role of multiple scattering and earth sphericity on 
the computation of photodissociation rates near dawn and 
dusk and subsequent effects on diurnal variations of strat- 
ospheric trace species. 

The document is organized into three technical sec- 
tions and six supporting appendices. Section 2 describes 
the results of 1-D model sensitivity studies of stratos- 
pheric ozone perturbation for a hypothetical fleet of high 
altitude aircraft. The implication of three models of OH 
chemistry and two choices of eddy diffusion profiles are 
discussed. In Section 3, approaches to formulation of 
multidimensional models are reviewed with emphasis on two- 
dimensional zonal-mean modeling. Difficulties associated 
with parameterization of eddy transport terms in Eulerian 
approaches are discussed and the Generalized Lagrangian 
Mean (GLM) formalism is presented as an alternative method- 
ology. Several crucial areas that require further investi- 
gation to practically apply GLM to zonal modeling are 
identified. In Section 4, possible approaches and 


approximations are examined for incorporating both spheri- 
city and multiple scattering in diurnal calculations for 
free radical species with emphasis on behavior near dawn 
and dusk. Diurnally dependent photodissociation rates and 
species concentrations are evaluated using a first order 
technique and an approach to incorporate a priori informa- 
tion on diurnal variations of photochemically active 
species within occultation based inversion algorithms ■ is 
discussed. A general summary is included in Section 5. 

We would like to thank R. Specht for his responsive 
programming support during the course of this work and 
K.K. Tung for fruitful discussions during the preparation 
of Section 3 of this report. 
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2 . 


MODEL SENSITIVITY STUDIES 


2 . 1 Overview 

Earlier model calculations (Crutzen, 1970; Johnston, 
1971; McElroy et al . , 1974; CIAP , 1974) indicate that a 
potential fleet of operations of supersonic aircraft as 
contemplated by the United States in 1970 (500 aircraft 
flying approximately 7 hours per day at 17-18 km) could 
lead to major reduction in the stratospheric column O 3 
abundance and thus cause an increase in the flux of ultra- 
violet radiation reaching the Earth's surface. This could 
result in a variety of environmental consequences, includ- 
ing a possible increase in the incidence of skin cancer 
(McDonald, 1971) . 

Removal of O 3 by aircraft injectant NO x radicals is 
primarily due to the pair of reactions 

NO + 0 3 N0 2 + 0 2 (2-1) 

followed by, 

N0 2 + O NO + 0 2 (2-2) 

The natural source for NO in the stratosphere is thought 

X 1 

to emanate from the reaction of 0( D) with N 2 0 (Nicolet and 
Vergison, 1971; Crutzen, 1971; McElroy and McConnell, 1971)/ 

N 2 0 + 0( 1 D) -> 2NO (2-3) 

Much of the revision in model predictions since 1976 
(Duewer et al., 1977; Turco et al., 1978) may be ascribed 
to a change in the rate constant for the reaction of NO 
with H0 2 / 
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(2-4) 


NO + H0 2 N0 2 + OH 

Recent measurements (Howard and Evenson, 1977) indicate 
that the rate constant for reaction (2-4) is faster than 
previously thought, by about a factor of 30. An increase 
in the rate constant for reaction (2-4) tends to shift the 
equilibrium in HO x from HO 2 toward OH. Below 30 km, 
removal of ozone by HO x radicals proceeds mainly by, 

OH + 0 3 -v H0 2 + 0 2 ( 2-5a) 

H0 2 + 0 3 + OH + 20 2 ( 2 -5b) 

We may note that reactions (2-5a) and (2-4) followed by 
photolysis of N0 2 , 

N0 2 + hv -> NO + O (2-6) 

do not affect odd oxygen removal. Thus addition of NO 
would reduce the catalytic role of HO x as described by 
equations (2-5a,b). Furthermore, an increase in OH causes 
an increase in the rate for the reaction, 

OH + N0 2 + M -> HN0 3 + M (2-7) 

with subsequent decrease in the concentration of the free 
nitrogen radicals (NO + N0 2 ) and reduction in the 
efficiency of the nitrogen cycle as a sink for odd oxygen. 

There is little doubt that OH plays a pivotal role in 
stratospheric chemistry and in the perturbed environment. 
While the kinetic data base for HO x reactions have been 
significantly improved over the past years, remaining un- 
certainties in H0 X chemistry still, perhaps, represent the 
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largest source of error for model predictions- In fact, 
a comparison between model calculations and observations 
reveals several significant discrepancies which might be 
attributed to errors in the calculated OH concentrations in 
the altitude region 15-35 km. It was argued elsewhere (Sze, 
1978; Sze and Ko (1980) that significantly lower stratos- 
pheric OH concentrations than those calculated by current 
models are needed to account for the observed gradients of 
CIO (Anderson et al . , 1979) and for the observed ratios of 
HN0 3 /N0 2 (NASA, 1977, 1979; McConnell and Evans, 1978) and 
HF/HC1 (Sze, 1978). 

Another area of major uncertainty concerns atmospheric 
transport of trace species. Current one-dimensional models 
parameterize vertical transport by the so-called eddy 
diffusion coefficients which were mainly derived from 
observation of N 2 O and CH 4 . These models therefore ignore 
horizontal transport, while the natural distribution of 
ozone exhibits significant latitudinal and seasonal varia- 
tions . 

In order to address the uncertainties associated with 
atmospheric transport and chemistry, a series of models will 
be investigated in an attempt to quantify the sensitivity 
of O 3 perturbations to different models of NO x injection by 
aircraft operations. Our approach will emphasize the 
uncertainties in the stratospheric OH distributions and 
their implications for perturbation studies. 

2 . 2 Model Results 

The results presented in this section are calculated 
by a one-dimensional model (Sze, 1978; Sze and Ko, 1979) 
with the rate data for oxygen-nitrogen-chlorine reactions 
taken from NASA (1979) , while those for HO x reactions are 
discussed in the text. The aircraft emission characteris- 
tics corresponding to different emission indices, fleet 
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size and engine types are summarized in Appendix F. The 
model atmosphere was taken from U.S. Standard Atmosphere 
Supplement (1966) . The diffusion coefficients (K z ) for 
most studies presented here taken from Wofsy (1976), 
although in some model studies, we also consider other 
choices of K z . 

We shall consider three models (A, B and C) of OH 
chemistry that could have important implications for O 3 
perturbations associated with aircraft operations. The key 
characteristics of models A, B and C are defined in Table 

2-1. 

Model A uses the rate constants for key H0 X reactions 
as recommended by NASA (1979) . Model B uses somewhat 


different 

rate constants 

for the following reactions. 


OH 

k 8 

+ H0 2 '* 

h 2 0 + 

°2 

( 2 - 8 ) 

OH 

k 9 

+HN0 3 

h 2 0 + 

N°3 

(2-9) 

ho 2 

k 10 

+ ho 2 

h 2 0 2 + 

°2 

( 2 - 10 ) 


k ll 





OH + O 3 -> H 2 O + O 2 (2-11) 

The rate constants kg, k g , k-^Q and k-j^ in model B are 
adjusted within their experimental uncertainties so as to 
give a lower OH concentration above 18 km. 

Model C is designed to investigate the possibility of 
reducing H0 X concentration, particularly in the lower 
stratosphere, by the reaction of OH with H 0 2 N 02 , 

k 12 

OH + H0 2 N0 2 H 2 0 + N0 2 + 0 2 (2-12) 


Reaction (2-12) could be a major sink for H0 X if the more 



TABLE 2-1 (a) 


* 

Reactions of Major Importance to Stratospheric OH Concentration 


REACTION MODEL A MODEL B MODEL C 


OH + HO 2 


h 2 0 + 0 2 

2.4(7) 

1.2(7) 

2.4(7) 

OH + HN0 3 

->■ 

H 2 0 + N0 3 

4.8(4) 

9.0(4) 

4.3(4) 

ho 2 + ho 2 


h 2 0 2 *f 0 2 

1.5(6) 

3.0(6) 

1.5(6) 

OH + 0 3 


h 2 0 + 0 2 

a < 940 

9.6(5) exp ( — — ) 

1.2(6) expC^jp) 

9.6(5) exp(-^-p) 

ho 2 + no 2 

M 

-»■ 

H0 2 N0 2 

0 

0 

NASA (1979) 

H0 2 N0 2 + hv 


ho 2 + N0 2 



Molina & Molina (1980 

H0 2 N0 2 

M 

-j- 

H0 2 + N0 2 



Graham et al. (1978) 

H0 2 N0 2 + OH 

-> 

H 2 0 + N0 2 + 0 2 

3.0(5) 

3.0(5) 

3.0(5) exp(-— ) 

H0 2 N0 2 + 0 


H 2 0 + N0 2 + 0 2 

6.0(2) 

6.0(2) 

6.0(2) 

ho 2 no 2 + Cl 

-»■ 

HC1 + N0 2 + 0 2 

1.8(5) 

1.8(5) 

1.8(5) 


Two body rates are in unit of m mol S 



00 


Reactions of Major Importance 


REACTION MODEL A 


OH + H0 2 

~T 

h 2 o t 0 2 

4 (-11) 

OH + HN0 3 


H 2 0 + N0 3 

8.0 (-14) 

ho 2 + ho 2 


H 2°2 + °2 

2.5 (-12) 

OH + 0 3 


h 2 o + o 2 

1.6 (-12) exp (-— -) 


M 



HO 2 + N0 2 

-*■ 

H0 2 N0 2 

0 

H0 2 N0 2 + hv 

->• 

ho 2 + no 2 



M 



H0 2 N0 2 


ho 2 + no 2 


H0 2 N0 2 + OH 


h 2 o + no 2 + o 2 

5.0 (-13) 

H0 2 N0 2 + o 


H 2 0 + N0 2 + 0 2 

1.0 (-15) 

H0 2 N0 2 + Cl 


HC1 + N0 2 + 0 2 

3 (-13) 


* 

Same as Table 2-1 (a) except that 

Two body rates are in units of cm^ molecule S 



3 (-13) 



recent cross-section data for H 0 2 NC >2 reported by Molina and 
Molinda (1980) are valid and if k 12 is faster than 1.2xl0 6 
m^ mol“l s -1 ( 2 x 10“12 cm 3 s -l) . 

Figure (2-1) shows the calculated OH profiles for 
models A, B and C. Note that the OH concentration in model 
B is about a factor of two smaller than that in model A 
above ^30 km, while the OH concentration is model C is 
about a factor of two to three smaller than that in model 
A below ^30 km. It should be noted that current measure- 
ments of OH are restricted above 30 km. Concentrations of 
OH below 30 km can only be indirectly inferred from other 
observed quantities such as the gradients of CIO and the 
HN 03 iN 02 and HF:HCl ratios. 

Figures (2-2) through (2-4) show the calculated pro- 
files of CIO, HNO 3 /NO;? and HF/HC1 , along with available 
observations. We may note that results from Model C pro- 
vide significantly better agreement with observations 
(Figures 2-2, 2-3, 2-4) than model A. 

2 . 3 Perturbation Studies 

For each model of HO x chemistry, we consider two 
different NO x injection altitudes, one at 15-16 km and the 
other at 20-21 km. Table 2-2 summarizes the calculated 
column perturbations which may result from the operation 
of 1000 aircraft, flying 7 hours per day at two different 
cruise altitudes. Figures (2-5a,b) and (2-6a,b) present 
the calculated local ozone and NO x perturbations associated 
with supersonic aircraft operation. 

Model A predicts increases of about 1.6 and 4 percent 
in column ozone for the 15 and 20 km injections respective- 
ly. On the other hand, model C predicts a small increase 
in column ozone for the low altitude injection but a fairly 
large decrease of about 6 percent for the high altitude 
injection case. The calculated column ozone perturbations 
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ALTITUDE (km) 


45 


40 


* 



CfO VOLUME MIXING RATIO 


Figure 2-2 

Calculated altitude profiles of CIO mixing ratio 

corresponding to 41° zenith angle. The observations 
from Anderson et al. (1979) and Menzies (1979) are 
included for comparison. 


ALTITUDE (km) 



0.1 1 10 100 
HNO3/NO2 CONCENTRATION RATIO 


Figure 2-3 

Calculated ratio of HNO 3 /NO 2 from models A and C. The 
computed ratios correspond to sunset (x = 90°) at 30°N. 
The data points are: 

• Evans et al. (1976) sunset 

o Harries et al . (1976) Harries (1978) noon 44°N 

□ Fontanella et al. (1975) sunset 45-50°N 
A Lowenstein et al . (1978) daytime 20-40°N 
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Figure 2-4 

Calculated HF:HC1 ratios corresponding to 30°N, equinox 
condition. The observations are: 

x-x Buijs et al . (1977) 

o-o Farmer and Raper (1977) 
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TABLE 2-2 


Sensitivity of Column Ozone Perturbation 
to Supersonic Aircraft Cruise Operation 


Number'. Correspond to A0 3 

15 km - 
injection 


in percent 

2 0 km 
in j ection 


Model A 


Model B 


+ 1.7 + 4.0 

+ .32 - 3.9 


Model C 


+ .25 


6.0 



ALTITUDE (km) 










ALTITUDE (km) 








from model B are similar to but somewhat smaller than those 
of model C. We may note that the range of calculated ozone 
perturbations span from +4 percent (model A) to -6 percent 
(model C) for the high altitude injection case. The range 
is considerably smaller in the low altitude injection case. 

It seems clear that sensitivity of column ozone 
perturbations to N0 X injection depends not only on the ' 
cruise altitudes but also on the subtle differences in H0 X 
chemistry as illustrated by models A, B and C. While model 
A uses the best estimates of rate constants for the H0 X 
reactions , it fails to account for several observations as 
discussed earlier. On the other hand, model C seems to 
give better agreement with observations (see Figures 2-2, 
2-3, 2-4) , although several of the rate constants for 
reactions involving HO 2 NO 2 need future laboratory studies. 
Until better kinetic data are available, it is difficult to 
rule out any of the models discussed here. 

The results presented above are based on Wofsy's (1976) 
eddy diffusion profile. Since ozone below 30 km is control- 
led mainly by dynamical transport, it is useful to find out 
how sensitive are model results to the choice of diffusion 
coefficients. Table 2-3 summarizes the results of models A 
and C which are calculated based on Chang's (1976) K z 
profile. The results seem to be not much different from 
those calculated by using Wofsy's (1976) K z profile. 

It should be recognized that our analysis on the sensi- 
tivity of column ozone perturbation to eddy diffusion 
profile is restricted to a one-dimensional model. A more 
realistic description of the ozone problem clearly requires 
a two-dimensional model, since the distribution of ozone as 
well as aircraft injection are essentially two-dimensional 
in nature. Thus the one-dimensional sensitivity study of 
K z profiles may be quite artificial, in the sense that the 
sensitivity of ozone perturbation to transport could be 


TABLE 2-3 


Sensitivity of Column Ozone Perturbations 
to Eddy Diffusion Coefficient 


Number Corresponds 


Model A 


Model B 


to AO^ (%) using 

15 km 
injection 

+ 1.3 

- 0.0 


Chang's (1976) K z 

2 0 km 
injection 

+ 1.5 


- 5.5 
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quite different in a 2-D model in which both horizontal 
and vertical transports are considered. 

2 . 4 Concluding Remarks 

We have performed a series of model calculations to 
study the sensitivity of column ozone perturbations to the 
injection of N0 X associated with the operation of super- 
sonic aircraft. Because of the coupling nature of hydrogen, 
nitrogen and chlorine chemistry, addition of N 0 X could 
either increase or decrease local stratospheric O 3 . A 
model with high background OH concentrations (e.g., model A) 
tends to predict an increase in column ozone for low and 
high altitude N0 X injections. On the other hand, a model 
with lower stratospheric OH (e.g., models B and C) tends to 
predict a small increase (vl%) in column ozone for the low 
altitude injection case but a fairly significant reduction 
in column ozone for the high altitude injection case. More 
reliable kinetic data are clearly needed to narrow the 
uncertainties discussed here. For instance, the rate 
constant for the reaction of OH with HO 2 NO 2 needs to be 
measured on a high priority basis. Accurate determination 
of this rate constant could either rule out or substantiate 
model C. 

Changes in local ozone and N0 X concentrations in the 
stratosphere may also affect the radiation budget. Recent 
calculations by Wang and Sze (1980) indicated that a 
doubling in NO x may perturb the stratospheric temperatures 
and surface temperatures by as much as +1 K and +.15 K 
respectively, mainly through redistribution of stratospher- 
ic ozone (Wang and Sze, 1980). While changes in stratos- 
pheric temperature by 1 K are unlikely to be important in 
stratospheric chemistry, changes in surface temperature by 
.15 K are considered to be quite significant when compared 
with surface temperature changes caused by other atmospheric 
trace gases (Wang et al. , 1976) . 
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3. TWO-DIMENSIONAL ZONAL-MEAN MODELING 


3 . 1 Background 

One of the tasks in atmospheric modeling is to attempt 
to simulate the behavior of a trace gas in the atmosphere. 
The local concentration of a trace gas is governed by the 
three-dimensional continuity equation 

{ Jt + Y- V)f = Q/p (3.1-1) 

where f(t,x) is the mixing ratio, v(t,x) the velocity wind 
fields describing the general circulation, p(t,x) the air 
number density and Q(t,x) is the local net production or 
loss (by chemical and/or physical transformation) of the 
trace gas. Equation (3.1-1) gives the time rate of change 
of f in the Eulerian description of fluid motion. The 
quantities f, v, p and Q are to be considered as Eulerian 
field quantities as functions of time and spatial location 
with coordinates x. 

In order to solve equation (3.1-1) for the specie 
concentration, one must be able to provide values of v and 
temperature T (T is necessary for calculation of reaction 
rates) as functions of space and time either by parameter- 
ization or by solving the system of dynamic equations. The 
atmospheric circulation is governed by the coupled system 
of dynamic and thermodynamic equations (cf. Lorentz , 1967) 
momentum equation 

dv -i 

-T- = -2ftxv =■ Vp - V$ (3.1-2) 

dt ~ ~ pM ^ 

thermodynamic equation 
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continuity equation 


(3.1-4) 


(3.1-5) 

d 9 

where -=-r- = ttt- + v *V is the total time derivative, 
dt 3t 

The above are to be considered as equations for the 
Eulerian field variables v, p, p and 0. The newly intro- 
duced symbols have the following meaning: 


ELG. = -nU* 

dt 


p V * v 


equation of state: ideal gas law 


p = pRT 


ft = angular velocity of Earth 
M = average mass of an air molecule 
p = pressure 

§ = the geopotential gz where g is acceleration due 
to gravity, z is geometrical altitude 
0 = potential temperature related to temperature 

T by 0 = T(^-) K with K = R/C^; R the gas constant 

and C the specific heat at constant pressure 
? . . T 

J = the diabatic influence with the heating rate 


per unit mass 
R = the gas constant 
T = temperature 

Note that we have left out frictional forces in the momen- 
tum equation ( 3 . 1 - 2 ) under the assumption that they are 
unimportant for large scale motion. 

The system of equations (3.1-2) to (3.1-5) is coupled 
to the specie equation through the J term which depends on 
distribution of gases such as 03 , CO 2 / N 2 O and CH 4 in the 
atmosphere. Thus, in principle, equations (3.1-1) through 
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(3.1-5) must be solved simultaneously as a system. 

In practice, the system of equations presents a 
formidable numerical problem and put enormous demand on 
both computer core memory and computation time. This is 
particularly true if one is interested in a realistic 
chemical scheme in order to simulate the distribution of 
the various species in the atmosphere. Besides, the set of 
exact equations also simulates phenomena of little interest 
for large scale motions. The following physical assump- 
tions are usually adopted: 

A) Replacement of the vertical momentum equation by 
the hydrostatic equilibrium condition, i.e., the 
pressure gradient force is balanced by geopotential 
term) 

= -M p g (3.1-6) 

It is observed that motions due to deviation away 
from hydrostatic equilibrium are restricted to 
oscillations about the equilibrium state with 
time scales of order hour. The adoption of 
equation (3.1-6) effectively filters out vertically 
travelling sound waves. 

B) Discard terms containing the vertical velocities 
in the remaining two components of the momentum 
equation. This approximation appears to be justi- 
fied because of the fact that the vertical compon- 
ent of the velocity is about two orders of magni- 
tude smaller than the horizontal components. 

C) Replace r in the resulting equation by a where r 
is the distance from the Earth’s center; a is the 
radius of the Earth. 

The above assumptions lead to a system of equations 
usually referred to as the primitive equations. Next, it 
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would be desirable to write the momentum equation in scaler 
form. For this purpose, it is convenient to introduce a 
coordinate system where the pressure p is used as the ver- 
tical coordinates in conjunction with A the longitude and 
4> the latitude. In this coordinate system, the components 
of the velocity vectors are 


u 

v 


= a cost}) 


dA 

dt 


a 


d<(> 

dt 


w 


dp 

dt 


The continuity equation takes the form 


3u 


a coscj) 3A 


, 1 3 / j\ - 3w 

+ t ( v cos* + — 

a cost)) 34> Y 3p 


= 0 


(3.1-7) 


This is sometimes written in the form 

V * v = 0 (3.1-8) 

However, equation (3.1-8) holds only in pressure coordin-* 
ates . Since equation (3.1-7) is derived using only the 
hydrostatic equation, it is more general than the incom- 
pressible fluid assumption. Finally, from equation 
(3.1-7), one can derive (cf. Lorentz, 1967) 


ds 

dt 


3s, 1 3 , ^ 

3t a cos4) 3 A [us) 


+ 


1 3 

a cos4> 3 4> 


(v cos4> 



(ws) 


(3.1-9) 


for any scaler s. 

In this coordinate system, the primitive equations are 
(cf. Lorentz, 1967; Holton, 1975) 
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3f . 1 

T 


3t a coscf) 3A 


8 (fu) + 1 


a coscf Xf (fv cos + af = “ 


(3.1-10) 


9u x ' 1 3 , 2 V , 1 

+ r* -7TT- ( U ) + 


St a coscj) 3A 


a cos*" sf (uv cos< t>> + s| < uw > 


uv tan(j> . , cr 

- — - 2[2v sm<j) - — -2- 


9 z 

a coscf) 3 A 


(3.1-11) 


= 0 


3 v 


+ 


9t a coscf) 3A 


3 (vu) + 1 


a coscf) 3cf) 


3 , 2 . » 

(v COS(j)) 


+ (vw) 


J tancj) , ^ , q 3Z 

+ + 2Qu sinrfi - — 7 — = 0 

a a 3cp 


(3.1-12) 


30 + 1 


3t a coscf) 3A 


' (0u) + 1 


a cos* f* (0V cos *) + sf < 6w > = J 


(3.1-13) 


3u 1 


a coscf) 3 A a coscf) 3cf> 


9 (v coscf)) + = 0 

0P 


(3.1-14) 


p = p R T 


(3.1-15) 


3 z 

g 3? 


1 

pM 


(3.1-16) 


In the above equations, the terms uv tancf) and u_ 2 tancj) are 

a a 

curvature terms arising from the non-Euclidean nature of 

the coordinate frames. The geopotential term - ~ and 
1 30 , . a 9 <J> 

-T co s (j, aX has been written out explicitly with $ = gz. 

In this arrangement, equations (3.1-10) through 
(3.1-13) are prognostic equations for f, u, v and 6 
respectively, while equations (3.1-14) and (3.1-15) can be 
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considered as diagnostic equations for w and p. Since p is 
being used as an independent variable, equation (3.1-16) 
serves the purpose of transforming between geometrical 
altitude and the pressure coordinate. 

We will now discuss ways of reducing the system of 
primitive equations to 2 dimensions for numerical modeling. 

Taking any of the dependent Eulerian variables in the 
primitive equations, e.g., f, by averaging over one of the 
spatial coordinates, one obtains f (t,y) where y is the 
vector representing the two remaining coordinates. 

The corresponding equations governing the averaged 
variables can be obtained by applying the averaging opera- 
tion to the primitive equations. In atmospheric modeling 
studies, one is interested in averaging over the longitud- 
inal coordinate to obtain f as a function of time, latitude 
and altitude, which can be viewed as the zonal-mean concen- 
tration. However, the exact physical interpretation of f 
must depend on the averaging process involved. 

In all of the 2-D zonal-mean models currently in use 
(cf. Hidalgo and Crutzen, 1977; Louis et al., 1974; Whitten 
et al . , 1977; Borucki et al. r 1976; Vupputuri , 1973; 

Harwood and Pyle, 1975) the Eulerian-mean 
operation is used exclusively. In the context of the 
present discussion, the most important feature in a 
Eulerian-zonal-mean model is the treatment of tracer 
transport. In these models, the flux into or out of a 
volume element is comprised of two terms: the flux due to 
advection by the zonal-mean velocity and the contribution 
due to the so-called "eddy fluxes". Several studies 
(Sawyer, 1965; Mahlman, 1975; Matsuno, 1979; Plumb, 1979) 
have found the above feature in the Eulerian-mean model to 
be undesirable and further work (Andrew and McIntyre, 1978; 
Dunkerton, 1978; Matsuno and Nakamura, 1979; Matsuno, 1980; 
Mahlman et al. , 1980) eventually led to the development of 


27 



the generalized Lagrangian-mean (GLM) theory. As demon- 
strated by Andrew and McIntyre (1978), the flux of a trace 
gas in the GLM formulation is given as an advection term by 
a GLM velocity. Further studies by Dunkerton (1978) showed 
that this GLM velocity could be approximated from consider- 
ation of atmospheric heat balance. 

It is the purpose of this section to discuss the 
possibility of applying the GLM theory to numerical zonal- 
mean modeling of atmospheric trace gases. In the next 
section, we will begin with a review of the current 
Eulerian-mean models in order to put the problem into 
proper perspective. 


3 . 2 Eulerian-Mean Model 

In the Eulerian-mean model, the averaging process is 
performed over Eulerian field quantities by integrating 
over one of the coordinates while holding all other coor- 
dinates fixed. In the case of the atmospheric models, the 
integration over the longitude is done along constant lati- 
tude circles, altitudes and time. Thus, given g ( t , c)> , p , X ) , 
(where t is time, <J) the latitude, p the pressure height 
coordinates and A the longitude) , one obtains 


g (t , <J> ,p) = - 

I dA 


(3.2-1) 


It is convenient to define g', the deviations from the 
mean , by 


g' (t ,<J> ,p. A) = g (t,<J>,p) - g (t,<f>,p,A) (3.2-2) 


It follows from equation (3.2.1) that 


g 1 (t ,<f> ,p. A) 


(3.2-3) 
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and 


gh = gh + g'h 


(3.2-4) 


When the Eulerian-mean operation is applied to 
equations (3.1-10) to (3.1-16), using the fact that the 
" — " operation commutes with the coordinate differential 
operators , we have 


||- + — x 4 (f V C0S(f>) + ^ (f w) = (%- F (3.2-5) 

9t a coscf) 3 Y 9p p' f 


9u . 1 9/— — ,2-'- 9/— — x 


+ 


9t a cos^cj) 9cj) 


27 k ( u v cos z <f>) + ^ (u w) - 2^v coscj) = -F 


u 


(3.2-6) 


9 v 1 
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9 -x 


9t a cost)) 9 4> 


(v cos 4> ) + ^ ( v w ) + 


-2 

u tan<j) 


(3.2-7) 


+ 2 u ft sin4> - — |x = -F 
Y a 9 4> v 


H + 5 XXX 4 (e v cos * ) + (e w) = J - F e (3 - 2 - 8) 


1 9 f— .. , 9w 

X XT' ( V COS(M + TT— = 0 

a cos4) 94> ■ 9p 


(3.2-9) 


H = RT (3.2-10) 

P 


9 z 1 RT 

g 9p Mp " Mp 


(3.2-11) 
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where F f = - gf (f'V cos.),) + g| (f’W) 


F = 
u 


a cos (J) 


27 3? (V ' U ' COS ^ + 3? (u ' W,) 


F = — t *4- (v'^cosij)) + (v'w 1 ) + — — n< ^ 

v a costj) 8cf) v 8p a 


F q = — 1 - - ~r tv ( 0 ' v ' cos 6 ) + (0'w') 

0 a COS(p d(p dp 


are the eddy flux divergence terms. Note that in equation 

(3.2-6) , we have incorporated the curvature term ( uv ^ an( ^ ) 

8 a 
into the term. It should be emphasized that the eddy 

flux terms are not constrained by the system of equations 

and must be specified or parameterized in terms of other 

variables . 

In zonal-mean models, a further physical assumption is 
usually adopted. From scale analysis, it can be shown (cf. 
Holton, 1975) that equation (3.2-7) reduces to 

2ft sin0u = SL (3.2-12) 

a o(p 

to lowest order of a Rossby number expansion. This is the 
geostrophic assumption which provides a good approximation 
for midlatitude dynamics though some authors (Harwood and 
Pyle, 1975) argued that the approximation is reasonable up 
to about 5°. The geostrophic assumption filters out occur- 
rance of gravity waves and other equatorial waves that 
cause the quasi-biennial oscillations in the stratosphere. 
Equation (3.2-12) is awkward to use as z is not one of the 
dependent variables. One can differentiate the equation 
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with respect to p and use the hydrostatic equation to 
obtain the thermal wind relation 



(5:) (E ) K 

2 52 sin<|> 


1 _30 
a £(J) 


(3.2-13) 


to replace equation (3.2-7) . This eliminates F v as an 
extra variable. However, one still needs ways of determin- 
ing F^, F^ and F Q to close the system of equations. 

All existing zonal-mean models use the method of para- 
meterization via eddy diffusion tensor originally proposed 
by Reed and German (1965) . Using mixing length type argu- 
ments, Reed and German argued that for a quantity x which 
is conserved along its flow, the eddy fluxes are related to 
the gradient zonal-mean quantity x via 


( 


v x 


W ' X 





(3.2-14) 


where K is a symmetric tensor. Following the procedure 
similar to the one suggested by Reed and German (1965) , 
Luther (1973) analyzed the heat transfer, temperature and 
wind variance data of Oort and Rasmussen (1971) and derived 
a set of K tensors as a function of time and space. It 
should be noted that the data only covered part of the 
northern hemisphere. Extrapolations were used to deduce 
the K's at places where there is no data based on results 
of Newell et al., (1966) and Wofsy and McElroy (1973) . 
Furthermore, the K's for the southern hemisphere are 
obtained by reflecting the northern hemispheric values in 
the appropriate seasons. The set of K's from Luther (1973) 
is by far the most complete with a set of values for each 
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month covering the meridional plane from the surface to 
60 km. In almost all of the zonal models, this set of K's 
provides the basis for eddy transport. We would like to 
emphasize here the limitation of this approach: 

(i) The original argument of Reed and German (1965) 
was presented to treat turbulence diffusive type 
motions. Thus, the approach may not be appropri- 
ate if there is organized wave-type motion in the 
zonal direction. This has been demonstrated in 
studies of 3-D circulation by Mahlman (1975) , 
Matsuno and Nakamura (1979) , Plumb (1979) and 
Matsuno (1980) in the case of planetary wave 
motions. To date, no satisfactory justification 
for applying the eddy diffusion theory to 
stratospheric motion has been presented. 

(ii) The argument only applies to conservative flow. 
Thus, strictly speaking, it can only be applied 
to inert tracers or to potential temperatures in 
adiabatic flow. However, in zonal-mean models, 
it is applied to chemical species with finite 
chemical lifetimes and to the potential tempera- 
ture where the condition of adiabatic flow is 
not strictly satisfied. 

The treatment of presents special problems and will 
be discussed later. 

We will now examine some of the zonal-mean models and 
discuss the various approaches adopted. The models can be 
separated into two groups depending on the number of the 
primitive equations it attempts to solve. In the first 
group of models, only equation (3.2-5) is solved as a prog- 
nostic equation. The variables necessary as input for the 
equation, i.e., v, w and T (necessary for calculation of 
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reaction rate constants) are parameterized. In addition, 
a chemical scheme has to be set up to calculate the produc- 
tion and loss term Q. The models of Hidalgo and Crutzen 
(1977), Louis et al . , (1974), Widhopf (1975), the NASA Ames 

Model (Whitten et al., 1977; Borucki et al . , 1976) and the 
Meteorology Office Model (Hinds, 1979) fall into this 
category . 

Strictly speaking, the v, w and T specified must satis- 
fy equations (3.2-7) through (3.2-9). The circulations 
used in most of these zonal models are based on the work of 
Louis et al. , (1974) who computed the circulation patterns 

for each season by solving the continuity equation (3.2-9) 
and energy equation (3.2-8) using the compiled observations 
of local meridional temperature distribution and heat trans- 
fer rates. However, interpolations and adjustments are 
usually made subsequently by the individual authors. In 
those cases, w is obtained by solving equation (3.2-9) once 
v has been chosen in order to assure mass conservation. In 
models that use pressure coordinates, the hydrostatic 
equilibrium assumption is implicit in equation (3.2-9) . In 
some models, e.g., the NASA Ames Model, geometrical alti- 
tude is used instead. In this particular model, the zonal- 
meaned continuity equation takes the form 


9_P , 1 9 

3t a cos6 3d 


(p v COS0) 



(p z) 



(3.2-15) 


where F = — — -r ttt- (p'v 1 cosQ) + ttv (p'z') and z = is 

p a cosd 3d K 3 t dt 

the vertical velocity. Equation (3.2-15) is then reduced 
to 


1 3 

a cosp 3d 


(p v co s 0 ) 


+ 



z) 


0 


by assuming that = 0, and F 
these additional assumptions are 


0. It is not clear how 
to be justified. 
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We will next discuss how to obtain Q/p once a photo- 
chemical scheme has been set up. The term Q/p consists of 
sums of terms each of which take one of the following forms 



p k ij (T) f i f j ; p 2 k ijk (T)f i f j f k 


for photolytic processes and two-body and three-body reac- 
tions respectively. Note that the temperature dependence 
of the reaction rates and k^^. are explicitly displayed. 

The longitudinal behavior of and f^ arise from the 
diurnal effect since each longitude is at a different local 
time at any instant in time. A similar problem arises in 
the 1-D model and is usually taken care of by putting 
diurnal averaged photolysis rates in the chemical scheme so 
that the resulting species concentrations represent diurnal- 
ly averaged values. In addition, the longitudinal 
distribution of the gas may be affected by wave motion in 
the zonal direction. No attempts are made to treat such 
effects. In the zonal averaging process, the term Q/p is 
obtained by taking the sum of the corresponding products of 
the zonal-meaned quantities. Eddy flux terms are ignored. 

Finally, it should be noted that because of the uncer- 
tainties in the input data, most authors adjust the wind 
fields and eddy diffusion tensors to obtain agreement be- 
tween calculated results and observations as a way to cali- 
brate their model. This must be taken into account in 
model validation studies. Also, as the input data are 
based on observations of the present atmosphere, such 
models are not suitable for large scale perturbation studies 
since there is no way to predict future circulation 
patterns . 

The next group of zonal models actually solve the sys- 
tem of equations simultaneously. The approach will be 
outlined below. For a more detailed treatment, the reader 
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is referred to the work of Harwood and Pyle (1975) and 
Vupputuri (1973) . 

Equations (3.2-6) and (3.2-8) can be written as 


9u 

9t 


A 


(3.2-16) 
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3t 


B 


(3.2-17) 


where A and B do not contain any time derivatives. If one 
differentiates equation (3.2-16) by p and equation (3.2-17) 
by 4) and takes the appropriate linear combination of the 
resulting equation, by virtue of the thermal wind relations 
equation (3.2-13), one obtains 


3A + (mXpo) 1 SB = 
20 , sincj) a 3cj) 


(3.2-18) 


Note that there is no time derivative in equation (3.2-18) . 
Equation (3.2-9) implies the existence of a function ip 
where 

- = _ _1 9 ± 

cos (}) 9p 

(3.2-19) 

- 1 3 ^» 

a cos4> 34> 

Using equation (3.2-19) to eliminate v and w in equation 
(3.2-18), we obtain a second order partial linear differen- 
tial equation for ip where the coefficients are functions of 
u, 0, J, , Fq and their spatial derivatives. The system 
of equations can now be solved as follows. Give f, 0, u, v 
and w at one instance in time, one can compute Q/p , J, F f , 

F„ and F and solve the prognostic equations (3.2-5), 
u u 

(3.2-6) and (3.2-8) for f, u and 0 at a later time. With 
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these new values, one solves equation (3.2-18) for \}j with 
appropriate boundary conditions and generate new v and w 
via equation (3.2-19). The terms Q/p are to be calculated 
as discussed before. The treatment of J is similar in that 
the term is obtained by using zonal-meaned values in the 
expression while ignoring eddy flux terms. and Fg are 

to be parameterized by eddy diffusion tensors. The treat- 
ment of F is more difficult because the u momentum is 
u 

definitely not conserved in the flow. Attempts to para- 
meterize F^ in terms of potential vorticity and potential 
temperature has little success (Green, 1970; Wiin-Neilsen 
and Sela, 1971) . Vupputuri (197 9) employed an ad hoc para- 
meterization while Harwood and Pyle (1975) deduced the 
momentum fluxes from Nimbus V SCR data. Until a more 
satisfactory theoretical approach is available, the satel- 
lite data approach seems to be preferrable at this stage. 

Thus, even in the case when one actually computes the 
circulation and temperature, one still has to depend on 
observations to parameterize the eddy flux terms. It is 
also noted in the calculation of Harwood and Pyle (1977) 
that the fluxes due to the mean circulation are often in 
opposite direction to eddy fluxes. Thus, the net transport 
is the difference between the two terms. This suggests the 
act of splitting the transport term into the advection and 
eddy fluxes components may not be appropriate. It is the 
desire to get away from eddy flux formulation that led to 
the development of the generalized Lagrangian-mean zonal 
models . 

3 . 3 Generalized Lagrangian-Mean (GLM) Zonal Models 

An alternative to the Eulerian approach of fluid motion 
is the classical Lagrangian approach. In such an approach, 
one is interested in how a quantity changes as one follows 
the motion of the fluid. More specifically, the property 
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g is given by the function g(t,s) where s(t,x) is the posi- 
tion vector of a fluid particle labelled by x. This method 
has been used to study the general circulation in three 
dimensions in several studies (cf. Riehl and Fultz, 1957; 
Krishnamurti , 1961; Danielson, 1968; Mahlman, 1969, 1973). 
Further developments in an attempt to apply the Lagrangian 
type dynamics treatment to zonal-mean quantities led to the 
development of the Generalized Lagrangian-Mean (GLM) theory 
(Andrew and McIntyre, 1978; Kida, 1977; Dunkerton , 1978; 
Matsuno and Nakamura, 1979; Matsuno, 1980). 

It should be noted that the GLM zonal approach is 
actually a combination of Eulerian-Lagrangian descriptions 
of fluid motion. One actually starts with the Eulerian 
equation (3.1.1) and applies the GLM. operator obtaining 

9f L v L 9f L _ -L P- + Q7p L (3.3-1) 

dt a 3(j) 

where f L (t,(j>,p) is the Eulerian field re~ultincr from the 
application of the GLM operator to f (t,<f>,p,A) ; v L and w L 
are the components of the GLM velocity fields obtained by 
averaging the Eulerian velocity field. Before discussing 
the definition of the GLM operator, it should be noted that 
the tracer transport in equation (3.3-1) is given by advec- 
tion only, thus, eliminating any possibility of cancellation 
between advection and eddy terms . 

The GLM operation can be defined by the following set 
of equations (Andrew and McIntyre, 1978). Given an Eulerian 
field g(t,x) we defined g^ 1 by 

g L = g(t,x + £ (t,x) (3.3-2a) 

where E, satisfies 
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(ft + y L .V) (x + g(t,x)) = v(t,x + §) 


or 


(ft + y L * v ) § (t»?) 


v(t,x + £) - V L 


(3.3-2b) 


and 


5 (t 


(3. 3-2c) 


In the above equations, " " denotes Eulerian averaging 

along the longitudinal coordinates, £(t,x) is the 
disturbance induced displacement field, v(t,x) is the fluid 
velocity field (an Eulerian field quantity) and v L is the 
GLM velocity field obtained by applying L to v L to v, i.e., 

v L = v(t,x + g(t,x)) . Note that from equation (3.3-2a), 

instead of averaging over constant latitudes and altitude, 
the averaging is done following the fluid particles to 
their displaced position. This gives a Lagrangian flavor 
to the approach though £(t,x) does not necessarily have to 
represent actual fluid motion. In the formulation, v L and 
E, are implicitly defined by equations (3.3-2b) and (3.3-2c) 
and can be obtained by solving the equations given v(t,x) 
together with appropriate initial conditions. And as long 
as equations (3.3-2b) and (3.3-2c) are satisfied, equation 
(3.3-1) can be obtained by applying ^ to equation (3.3-1) . 

To date, most of the work in GLM analysis has been 
done from a three-dimensional point of view, i.e., the 
motion (usually wave type) in 3 dimensions is specified, 
the properties of wave transport are then studied in the 
context of GLM theory (cf. Matsuno and Nakamura, 1979; 
Matsuno, 1980, Mahlman et al., 1980). 
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In a zonal model, one would like to integrate equation 

(3.3-1) in time to obtain £ as a function of time. Thus, 

■™L — L 

it is necessary to obtain v and w without a priori 

knowing the 3-dimensional general circulation. Secondly, 

one has to be more specific in the definition of E, so that 

a physical interpretation can be given to the GLM operation 
-L 

and that f can be compared to observations. We will 
discuss these problems in what follows. 

(a) Interpretation of E, and L 

In the formulation of Andrew and McIntyre (1978) , the 
definition of £ is purposely left vague in order that the 
formalism can be applied to any averaging operations. The 
quantity £(t,x) is loosely referred to as the disturbance 
induced displacement vector. One can imagine the situation 
in which a parcel of air molecules is initially at rest on 

o 

a latitudinal circle at constant altitude [x] . (Note that 

O ~ 

[x] stands for a latitude circle at constant altitude) . At 
t = t c , the molecules start executing motions given by 

O 

£(t,x), i.e., the coordinate of the particle initially 

^ o o o 

located at x is given at time t by x + £( t,x) . Then 
according to equation (3.3-2a) , instead of averaging over 

O 

the particles that happen to be located at x at time t, the 
averaging is done following the fluid particles as they 

o o 

travel along the trajectory £(t,x). If £(t,x) satisfies 
equations (3.3-2b) and (3.3-2c), then the averaging pro- 
ceeding qualifies as GLM operation, and v^ can be inter- 
preted as the motion associated with the center of mass of 

o 

the air parcel originally located at x. This interpreta- 
tion has been presented by Andrew and McIntyre (1978) using 
their mechanical analog and is essentially the same as 
interpretations presented by other authors (Dunkerton, 

1978; Mahlman et al., 1980). 
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The above interpretation of v L is contingent upon the 

o 

fact that £ (t ,x) satisfies equations (3.3-2b) and (3.3-2c). 
In general, this will not hold for Lagrangian trajectory of 

o 

the fluid. However, since £(t,x) is defined by the current 

~ ~ o 

positions of the particles as well as x, one can redefine 

£ as £(t,x c ) so that equations (3.3-*2b) and (3.3-2c) are 

satisfied (Tung, private communication) . The details of 

this manipulation will be presented in Appendix E. 

-L 

Under the initialization process, f can be inter- 
preted as the average mixing ratio of the air molecules 
that would have been located at x D in the absence of the 
disturbance induced displacement. Similarly, Q/p L is the 
average production rate of the same set of molecules and 
v L and w L are the components of the velocity associated 
with the center of mass of that air parcel. 

Note that a parcel of air has a tendency to spread out 
because of diffusion. Thus, the magnitude of £(t,x 0 ) would 

get progressively larger. It will then become increasingly 

-L . . 

difficult to relate f to actual observations as the indiv- 
idual molecules that enter into the averaging process are 
distributed over the entire atmosphere. The same applies 
to the computation of Q/p L in the sense that Q/p L at some 

particular latitude and height may differ significantly 

— L 

from the corresponding products of the f s at the same 
point . 

(b) Determination of v L and w L 

In order to use equation (3.3-1) as a prognostic 
equation, one needs the values v^ and w . If the three- 
dimensional motion is known, one can conceivably solve 
equation (3.3-2b) and (3.3-2c) for E, and v L with choice of 
initial conditions. However, if one is interested in a 
zonal model than can include feedback mechanisms between 
chemistry and dynamics, one would have to be able to obtain 
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v L and w L without a priori knowledge of the motion. 

One possible approach is to try to obtain v L and w L 
from the GLM momentum equation. Unfortunately, as demon- 
strated in Andrew and McIntyre (1978) , the price paid for 
the simplification in equation (3.3-1) is that L no longer 
commutes with the coordinate differential operators. For 
example, the GLM flow of an incompressible fluid satisfying 
V.v = 0 is not non-divergent , i.e. V.v L ^ 0 , in contrast to 
the case of Eulerian mean flow. Another result of the non- 
commutability of the operator L and V is that the GLM 
momentum equation contains terms involving i.e. 

Vp L = Vp L + (term involving E, and derivative of ^) 


These terms are reminiscent of the eddy flux terms in the 
Eulerian mean formulation and require parameterization. 

Dunkerton (1978) suggested an alternative approach to 
the problem from analysis of the thermodynamic equation. 
Using a linearized version of GLM theory for steady waves, 
Dunkerton showed that under certain conditions 

-L Q 

w = r d - r (3.3-3) 


1 _9 

a coscf) 9(J> 



cos<}) ) 


9w L 


0 


(3 . 3-4) 


where Q is the Eulerian mean diabatic heating rate, 1 ^- T 
is the lapse rate. The 3-D equation for potential tempera- 
ture can be written 

(|^ + v.V)T = Q (3.3-5) 


where Q is the diabatic heating. Upon GLM averaging one 
obtains 
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(3.3-6) 


9T L v L 9T L -L 9T L 

at a 3 cf) W 9p 


o L 


Equation (3.3-3) can be derived from equation (3.3-6) if 
the following assumptions are made (Dunkerton, 1978) 


(i) 

9T L 

at 

negligible 

(ii) 

v L 9T L 

a 3<j) 

negligible 

( iii) 

-L ^ — 

Q ^ Q 


(iv) 

-L ^ — 
T o, T 

. . . st l „ 

s° that ^ a. 


Equation (3.3-4) is the continuity equation obtainable from 

2 

the GLM equation by discarding terms of 0 (£ ) . Using the 
heating rate of Murgatroyd and Singleton (1961) , Dunkerton 
solved equations (3.3-3) and (3.3-4) for v L and w L and 
obtained results that are capable of explaining the gross 
features of tracer transport. 

Although the assumption used in deriving equations 
(3.3-3) and (3.3-4) may not be valid in the whole atmos- 
phere, the result of Dunkerton is beautiful in its concep- 
tual simplicity in relating the circulation to diabatic 
heat balance. However, the application of the results to 
zonal modeling may not be straight forward. The situation 
is particularly difficult in dealing with finite amplitude 
motions . 

3 . 4 Concluding Remarks 

In this section, we have presented a brief review of 
the generalized Lagrangian mean description of atmospheric 
circulation. Most of the studies done to date in GLM 
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theory have been concerned with applying the technique in 
analyzing the transport due to specific wave motions [in 3 
dimension] in the atmosphere. In this section we have 
taken the view of trying to use equation (3.3-1) as the 
prognostic equation in zonal modeling. 

The zonal mean tracer equation (3.3-1) is attractive 
in that the only transport of tracer is due to advection by 
the Lagrangian mean velocity only. From our studies, we 
identified several crucial areas that need to be investi- 
gated in order that the GLM tracer equation can be profit- 
ably applied to zonal modeling. The first of these has to 
do with the basic interpretation of equation (3.3-1). As 
demonstrated in Section (3.3) and Appendix E, the distur- 
bance induced displacement vector can be related to the 
actual fluid motion. Then f^ can be interpreted as the 

average mixing ratio of the air molecules in an air parcel 

° ° 

originally located at [x] where [x] stands for the coordin- 
ates of the set of points on a latitude circle at constant 
altitude. Though the idea of keeping track of the identity 
of the molecules may be very useful in the analysis of mass 
transport, it may introduce complications in modeling of 
trace gases. As an air parcel is spread out thinner and 
thinner over an ever increasing region of the atmosphere, 
it becomes almost impossible to relate f L to the local 
observed concentrations. This suggests that one may want 
to relax the rigid relation between E, and the trajectory 
introduced in Appendix E and employ different initializa- 
tions during the time integration of equation (3.3-1) using 
each initialization for only short time intervals. 

An ideal zonal model should be able to compute the 
circulation given the ambient concentrations of the trace 
gas and the state of the atmosphere. Thus, it would be 
nice if there is some way to solve for v L and w L as 
functions of f L s and temperature. Such an approach will 
require research into the reformulation of the momentum 
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equation. By defining some kind of generalized differen- 
tial operators, it may be possible to cast the momentum 
equation in a form that has no explicit dependence in 

The above areas involve long term developments. On 
the more immediate term, one can still explore the merits 
of the GLM model by using v L and w L derived from given 
three-dimensional motion fields. The initial analysis 
could be applied to studies of inert tracers, thereby 
avoiding the question about Q/p^. Numerical experiments 
dealing with inert tracer release can be done. Such exper- 
iments will be more meaningful if the results can be 
compared with results from similar experiments using 
Eulerian mean models. One can then compare the center of 
mass motion obtained from the Eulerian model to see if it 
bears any resemblence to v^ and w^. 

One could also concentrate on cases with small distur- 
bances and make use of the results of Dunkerton (1978) to 
generate the wind fields. The heating rate of Murgatroyd 
and Singleton (1961) was used in the original analysis by 
Dunkerton (1978) . We feel that a more careful treatment of 
the diabatic heating could be used to generate v L and w^ 1 so 
that they can be used in subsequent tracer studies. Again, 
the study of such General Lagrangian Mean transport can be 
carried out most profitably in conjunction with parallel 
Eulerian mean model studies. One can take the results of 
Eulerian mean cell models which contain calculated v, w, 
and the temperature distribution. The temperature struc- 
ture can be used to generate the Lagrangian velocity field 
v^ and w L via Dunkerton * s approach. At the same time, one 
can use the Eulerian model to perform tracer experiments 

and try to determine if there is any relationship between 

— L — L 

the center of mass motion of the tracer and the v and w . 
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4. REMOTE SENSING 


4 . 1 Background 

Spectroscopic measurements using remote sensing tech- 
niques can provide information on the altitude distribu- 
of constituents in the upper atmosphere. Using a solar 
occultation approach, for example, inversion methods have 
been demonstrated to obtain ozone, aerosol, and nitrogen 
dioxide profiles from solar extinction data in the .38-1.0 
pm wavelength region (Chu and McCormick, 1979) . Since 
occultation measurements are made at definite times of the 
day, i.e., twilight or dawn for solar occultation, it is 
clear that diurnal variation of species must be considered 
in the interpretation of results. Particularly rapid vari- 
ation of photochemically active species (such as OH, H0 2 , 
CIO, NO, etc.) (see Logan et al. 1978) may introduce 
ambiguities in the treatment of species abundances deter- 
mined along a particular line of sight since such experi- 
ments are essentially looking through different local times 
and altitudes (Herman, 1979). Conversely, inversion 
results may provide a stringent test for diurnal model 
calculations . 

Additional constraints relate to the role of molecular 
multiple scattering on both transmitted radiances and 
dissociative fluxes required for photochemical modeling. 

As pointed out by Callis (1974) , molecular scattering of 
solar radiation can significantly affect photolysis rates 
in the lower atmosphere (z < 40 km) . A number of studies 
(Callis et al . 1975; Luther and Gelinas, 1976? Luther et 
al . 1978; Anderson and Meier, 1979) have examined multiple 
scattering effects on photolysis rates and species concen- 
trations employing the plane-parellel approximation of 
radiative transfer. Near dawn and dusk, however, the 
inherent curvature of the atmosphere becomes an important 
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consideration and realistic treatments of the effect of 
multiple scattering require radiative transfer calculations 
performed within the pertinent "spherical cap" geometry. 

In this section possible approaches and approximations 
are described for incorporating both sphericity and 
multiple scattering in diurnal calculations for free radi- 
cal species (e.g., CIO, ClNO ^ r HC^ / etc.) with emphasis 

on the behavior of concentrations at dawn and dusk. The 
effect of molecular scattering and surface reflection is 
included in the calculation of relevant photodissociation 
rates. In support of these aims, the following tasks have 
been performed: 

. investigation of optimal methods to evaluate the 
primary source function in the spherical shell 
geometry 

. formulation of the multiple scattering problem 
including surface reflection 

. evaluation of profiles of atmospheric optical 
properties at dissociative wavelengths 

. evaluation of single scattered dissociative fluxes 
and comparison with both pure absorption and 
multiple scattering results 

. evaluation of single scattered photodissociation 
rates and comparison with both pure absorption and 
multiple scattering results 

evaluation of single scattered photodissociation 
rates for the spherical shell atmosphere near dawn/ 
dusk 

. calculation of diurnal variations of radical 
species using single scatter photodissociation 
rates for the spherical geometry and comparison 
with plane-parallel pure absorption results 

. formulation of an approach to include higher orders 
of scattering and investigation of its convergence 
to multiple scattering results 

. formulation of an approach to incorporate a priori 
information on diurnal variations of photochemical- 
ly active absorbers within occultation based 
inversion algorithms. 
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'These results may be used to classify the atmosphere 
into regions where: (a) molecular scattering may be ignor- 

ed, (b) single scattering can adequately account for the 
effect of molecular scattering, and (c) higher order 
scattering must be considered. Such information is crucial 
to the development of computationally efficient model 
algorithms . 


4 . 2 Role of Earth Curvature and Scattering in Diurnal 
Photochemical Modeling 


Photodissociative processes play a significant role in 
current descriptions of atmospheric chemistry. Even the 
basic four reaction, oxygen-only model of Chapman (1931) 
includes two key photolytic reactions : 

0 2 + hv + 20 (4-1) 

0^ + h’O ->• 0 2 + 0 (4 — 2) 


The photodissociation rate, j (s ^) , of such processes is 
given by: 


j 1 ( z / 9 o ) 



K F x (z ' 0 


dX 


(4-3) 


where 



is the wavelength dependent photodissociation 
cross-section for species i (cm^) 

is the corresponding wavelength dependent quantum 
yield (nd) 

is the available dissociative flux at wavelength 
X, altitude z, and solar zenith angle 6 0 (photons 
cm ^ s 1 a ± ) . 


The time of day, t, determines the value of the local solar 
zenith angle 0 O as described in Appendix A. The dissocia- 
tive flux is given by the angle integrated specific 
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(4-4) 


intensity 1^ (photons cm 


“2 a~ 1 A - 1 


sr 1) 


F^(z,e„) = f i x (z,0 o ) an 


Q 


where 1^ is the solution to the general radiative transfer 
equation (Lenoble, 1977) : 


(ft.V)I A (r,ft) 


-3^ (r) [I A (r , ft ) 


J A (r,ft) ] 


(4-5) 


where r is a position vector, ft is a unit vector in the 
direction along which variations in the intensity are 
sought, 3^ (m ^) is the total extinction coefficient 
describing the loss in intensity due to both scattering and 
absorption processes, and is the source function deter- 
mined by scattered contributions to the intensity field. 

Solution of equation (4-5) for a given model atmos- 
phere requires an understanding of both: (a) the geometry 

of the atmosphere and (b) the relevant transfer processes 
(i.e., absorption, scattering) involved at a given wave- 
length. For example, in the simplest case, assuming a 
plane-parallel (PP) atmosphere, (appropriate for 0 O < 80°) 
ft reduces to the local normal n and dependence on r to that 
on altitude, z (km) . If additionally scattering is neglec- 
ted (i.e., pure absorption (PA)], the source function J, is 

d ^ 

identically zero and a simple form for F^ emerges: 


d , PP , PA 

F x 


F a T a PP (z,6 0 ) 


(4-6) 


where ttF, is the unattenuated solar flux at the top of the 
^ PP 

atmosphere and (z,0 o ) is the transmission function from 
the top of the atmosphere to altitude z along a path at 
solar zenith angle 0 O given by: 
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(4-7) 


T"(z,0 o ) = 

K . 

exp|-sec0 o Z [a^N^(z) ]j> 

i-1 

for K total absorbers where N^(z) is the vertical column 
density (cm - ^) f or absorber i from level z to space. The 
quantity in brackets is the corresponding vertical (or 
normal) optical depth, 

When scattering is included in the plane-parallel (PP) 
geometry, a variety of techniques are available to solve 
the radiative transfer equation (see Hansen and Travis, 
1974; Lenoble, 1977). However, near dawn and dusk, the 
large solar zenith angles characterizing these situations 
(0 O > 80°) necessitate consideration of the Earth's 
sphericity in evaluating path extinction of incident solar 
radiation. Thus, even in the pure absorption (PA) case 
described above, the transmission function (4-7) must be 
modified to account for Earth curvature. Appropriate 
modifications are discussed in §4.3. 

When scattering processes are additionally to be 
included, appropriate solutions to the general radiative 
transfer equation (4-5) must be sought. These will be 
described in §4.4. 

4 . 3 Optical Paths in a Spherical Shell Atmosphere 

4.3.1 Geometric Considerations 

As the solar zenith angle, 0 O , approaches 90° (i.e., 

near dawn and dusk) , equation (4-7) suggests that the 
transmission of all incident solar radiation (regardless of 
wavelength) approaches zero vary rapidly. This is contrary 
to experience. It is generally recognized that application 
of the plane-parallel approximation becomes inappropriate 
in these special cases and furthermore other physical 
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phenomena such as refraction enter the problem. To facili- 
tate evaluation of diurnally dependent dissociative fluxes, 
it is desirable to consider formulations which allow for a 
smooth transition between plane-parallel situations and 
those constrained by the spherical shell nature of the 
atmosphere. Ideally values of the solar zenith angle 
greater than 90° should be admitted to account for illumin- 
ation of higher altitudes (above the terminator) at predawn 
and post-dusk times. In this section, generalizations of 
the direct solar beam transmission function analogous to 
equation (4-7) are discussed. The significance of this 
function twofold: (1) when scattering is neglected, it 

determines the locally available dissociative flux, and (2) 
when scattering is treated, it additionally enters into the 
calculation of the primary source function for scattering. 


4.3.2 Air Mass Factor Formulation 


Transmission of incident solar radiation of wavelength 
A along a ray path s at solar zenith angle 0 O to the local 
normal n at altitude z will be given by the transmission 
function T: 


T a (z,0 o ) = exp [-t a (z, 0 o )] 


(4-8) 


where the slant path optical depth t a (z, 6 0 ) is defined as: 


T ^ ( Z , 0 o ) 


K 

V 1 

x °x 

1 



(z,0 o ) • 

n 1 (s' )ds’ 


(4-9) 


i th 

where: n x (s) is the number density of the i species along 

geometric ray path s (cm ) 

and ds ' is the incremental path length along the ray 
path from source to level specified. 

This expression allows for extinction due to K optically 
active species. The quantity in brackets may be defined 
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as the slant path column density (cm 2) Q f species i or: 


• (s(z,0 o ) . 

N 1 (z,0 o ) = J n 1 (s') ds ' (4-10 ) 


In the plane-parallel (PP) limit, the incremental path 
length, ds ' , and its projection on the local normal vector 
n, dz ' , are related by: 

ds 1 = sec 0 0 dz 1 (4-11) 

yielding : 

K i i 

t^(z, 0 o ) = sec 0 0 Z ct a JSH(z) (4-12) 

i 

where 

/ °° . 

n 1 (z ' )dz ' (4-13) 

z 

and 

^(z.e.) = sec e o N*(z) (4-14) 


Here N^(z) is the vertical column density of species i from 
level z to the top of the atmosphere along the normal 
vector. The ratio of slant path to vertical column density 
defines a plane-parallel slant path air mass factor for all 
i and z : 


M ( 0 0 ) = 
PP 


N 1 (z,0 o ) _ 


N 1 (z) 
v 


- = sec 0 


(4-15) 


Thus, from equation (4-12), the slant path optical depth 
may be determined from: 


(z , 0 0 ) 


I\ 

Z tV(z, 0 o ) 


(4-16) 
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where is the wavelength dependent vertical optical depth 
above level z given by: 

T*(z) = N*(z> (4-17) 

The approach employed in equation (4-16) may be applied to 
situations where the plane-parallel approximation is not 
applicable by formulating appropriate generalized func- 
tions , M 1 ( z , 0 o ) . 

4.3.3 The Chapman Function 

For grazing incidence (values of 0 O approaching and 
greater than 90°), the geometry of the problem is illus- 
trated schematically in Figure 4-1. The atmosphere is 
assumed to be spherically symmetric, that is the number 
density profiles are of the form: 

n^ = n^(r) (4-18) 

where 

r = R + z (4-19) 

and R is the planetary radius. (The value of R is assumed 
to be 6371 km, the Earth's radius at about 35°N.) Neglect- 
ing refraction, from simple geometric considerations: 

(R + z) sin0 o = (R + z 1 ) sin6 ' (4-20) 

and 

ds' = dz' sec8 ' = dz 1 [1 - (|^|, ) 2 sin 2 8 0 ] "^ (4-21) 

The desired slant column density at z along the path at 
solar zenith angle 0 O will be [from (4-10)]: 

N 1 ( z , 0 o ) = J n 1 (z ' ) [1 - (|±|,) 2 sin 2 e o r 2 dz' (4-22) 


52 



Assuming a number density distribution with constant scale 
height H 1 : 

n^(z') = n^ (z) exp (-z */H^) (4-23) 

and using 6' as the integration variable, (4-22) may be 
expressed as: 

N 1 (z,0 o ) = 


n 1 (z)H 1 X 1 sin0 


/ 


csc 2 0 ' expX 1 (1 - 


sin0 

sin0 


-)'dQ' 


= n 1 (z)H 1 Ch(X 1 ,0 o ) 

= n 1 (zjH 1 !^ (z,0 o ) (4-24) 

where Ch(X 1 ,0 o ) is the Chapman (1931) function and the 
parameter X 1 is : 

X 1 = (R + Zi/H 1 (4-25) 

This expression has been tabulated by Wilkes (1954) and is 

valid for 0 O £ 90°. Figure 4-2 compares the Chapman 

function (0 o ) with the plane-parallel air mass factor 
bp 

(i.e., sec0 o ) M pp (0 o ) for a uniformly mixed gas (02 or M) . 

(0 o ) assumes a constant scale height equal to 8.0 km. 
bp 

For 0 O > 90°, the Chapman function has been analytically 
extended using (see Figure 4-1) : 


N 1 (z,0 o >9O°, ® *> © ) = 2N 1 ( z o ,0 o = 90 °, © - ® ) 

-N 1 (z , 18O-0 o , © ) 


(4-26) 


although this approach is not appropriate as discussed by 
Swider (1964) . Note that Mp p and Mg p diverge as 0 O 
approaches 90°. From equation (4-16) the slant path 
optical depth may be expressed as : 
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Figure 4-2 

Comparison of air mass factors M(0o) for: (1) Chapman func 
tion, (2) Generalized Chapman function, (3) Numerical, (4) 
Numerical with Refraction, and (5) Plane Parallel. 





(z,e Q ) 


(4-27) 


= r ch ( x 1 , 0 „ ) 

i 

with appropriate choice of X 1 . Extensions to the Chapman 
function include work by Swider (1964) who investigated 
atmospheres with scale height gradients and Green and 
Martin (1966) who generalized the analysis for various 
density distributions. 


4.3.4 Other Analytical Treatments 

Swider (1964) noted that in the terrestrial atmosphere 
scale height variations with altitude occur and values of 
the Chapman function are most sensitive to parameters 
(scale heights and number densities) near the lowest 
altitude level encountered along the appropriate ray path. 
For zenith angles greater than 90°, this level is near the 
grazing altitude (or tangent height) z Q (see Figure 4-1) . 

By incorporating these considerations simplified approxi- 
mations have been derived by numerous authors including 
those by Smith and Smith (1972) from the work of Swider 
(1964) and Fitzmaurice (1964). For 0 O < 90°: 


M s P (z ' e 


< 90°) = 


(- 




i x 2 


X z ) exp (y ) erf c y‘ 


(4-28) 


Noting that Mg^(z / 0 o = 90°) (^r X^) 2 relation equation 

(4-27) can be used to derive an expression for 0 O > 90°: 


M s P (z ' e 


> 90°) = (£ x;) 


[2ni - n^exp(y 1 ) 

H 1 n 1 Z 

z z 


(4-29) 


erfc y 1 ] 


where : 


i 

n o 


number density at tangent height z Q 
number density at level z 
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hJ = scale height at tangent height 
= scale height at level z 

xl = (R+ Zo )/hJ 
X* = (R+z)/H^ 

y 1 = (X^/2) 2 jcos6 0 | 


and erfc = complimentary error function (=l-erf) . 
Expression (4-29) accounts for scale height variability in 
a realistic terrestrial atmosphere. Figure 4-2 compares 
air mass factors m| as a function of solar zenith angle 
with the constant scale height Chapman function. was 

bp 

evaluated using the scale height and total number density 
profiles of the U.S. Standard Atmosphere (1976). 


4.3.5 Numerical Forms (including refractive effects) 

Slant path column density N 1 (z f 0 o ) [and hence air mass 
M(z,0 o )] may be explicitly evaluated by numerical integra- 
tion of the line integral in equation (4-10) along the 
appropriate ray path. Advantages offered by this approach 
include capabilities to: (a) incorporate realistic number 

density profiles, n 1 (s') , for each species, and (b) include 
the effects of refraction within the atmosphere. 

Arbitrary number density profiles may be incorporated 
into the analysis by integrating equation (4-22) with 
suitable quadrature for a specified profile n 1 (z') taken as 
a series of concentric locally homogeneous shells. The 
resultant air mass factor M |p( Q o) using the number density 
profile cited above is also illustrated in Figure 4-2. For 
comparison, previously discussed air mass factors are also 
provided . 

The effect of refraction on optical path for zenith 
angles, 0 O > 90°, is schematically illustrated in Figure 
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4-3. In the absence of refraction, the appropriate optical 
path integral is along (b) at tangent height z 0 , for rays 
reaching altitude z at astronomical zenith angle, 0 O . The 
unrefracted tangent height is simply determined geometri- 
cally : 

z 0 = (R + z) sin ( 180-0 0 ) - R (4-30) 

for 0 O > 90°. With refraction, the apparent angle of 
arrival is 0£ and the minimum height of the refracted ray 
is Zq* For monotonically decreasing density profiles the 
minimum height of the refracted path is always greater than 
that for the unrefracted path, Snider (1975) : 

z' 0 - z 0 > 0 (4-31) 

Sample calculations of the error in the minimum height due 
to refraction equation (4-31) performed using the procedure 
described in Appendix B are shown in Figure 4-4 for values 
of z = 25, 35, 45, 55, 65 km. Note that the error is 
greatest for high altitudes at large zenith angles (low 
minimum heights). Since the column density for 0 O > 90° is 
most sensitive to the relevant number density near the 
tangent altitude [n Q , see equation (4-29)], the effect of 
refraction is to decrease the air mass. Neglecting refrac- 
tion, therefore, causes air mass to be overestimated 
(Snider, 1975) . 

4 

Calculations of refracted air mass factor M c (0 O ) 

bp 

follow a procedure based on Selby and McClatchey (1972) . 
Essentially, equation (4-22) is modified by a factor 
dependent on the phase refractive index of air, m(z) with 
the resultant air mass factor given by (see Appendix B) : 
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(Zo-z 0 ) ERROR IN MINIMUM HEIGHT (km) 


6.0 



Figure 4-4 

Calculated Error in Tangent Height Due to Refraction at 
Indicated Altitude. 
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1 


M 


Sp 


lz ' 9 - 1 

*} 


( 4-32 ) 


sin^G n \ ^dz' 


Results based on the procedure detailed in Appendix B are 
also illustrated in Figure 4-2. Note that the refracted 
air mass is less than that for a realistic profile ignoring 

3 

refraction (i.e. , M ) . This error for selected angles is 

bp 

given in Table 4-1. The error increases markedly for large 
solar zenith angles. 


4.3.6 Summary 

The differences among the various spherical approaches 

illustrated in Figure 4-2 is not as striking as that 

between the plane-parallel and spherical treatments. To 

avoid time consuming line integrals and yet retain the 

scale height dependence of the true atmosphere we have 
2 

retained the M ( 0 O ) approach (which is intermediate in 
bp 

value) . As applied to model calculations of solar beam 
extinction, optical paths are evaluated using equation 
(4-16) such that: 


t x (z ,e Q ) 

£ T 1 M 1 
i=l v 

(z, 0 o ) 


(4-33) 

where M 1 (z, 0 o ) = 

sec 0 o 

V 

0 

CD 

0 

0 

00 


II 
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CD 

V 

N 

s 

M? p (e 0 ) 

0o 1 

00 

0 

0 



In the current treatment three species (i.e., sources 
of extinction) are handled (i.e., K=3) : O 2 / 0 3 , and 
Rayleigh scattering (M) . Appropriate choices of scale 
height and number density profiles are used in (4-28, 4-29) 
depending on whether air mass factors for ( 0 3 * M) or 03 are 
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TABLE 4-1 


Computed Air Masses and Error in Neglecting Refraction 


' 6o 

M sp (e °) 

M sp (6 °> 

A ( % ) 

92.3 

126.4 

122.0 

3.7 

94.0 

700.5 

580.8 

20.6 

95.2 

4177.6 

2999. 

39.2 

96.4 

22030. 

15328. 

43.7 
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being calculated. Thus, the transmission function in a 
spherical shell atmosphere T^ p (z,0 o ) [analogous to equation 
(4-7)] becomes: 

Q - K . 

T^ p (z,0 o ) = exp -[Z T r" M (z,0 o ) ] (4-34) 

A i=i v 

and the dissociative flux due to attenuation (including 
scattering of the direct solar bean (but excluding multiply 
scattered contributions) is : 

F x' SP = F X T x P(z ' 0 » ) (4-35) 

For ozone, a similar approach is adopted (valid above 
25 km) modifying the ozone scale height to be equal to 2/3 
the ambient pressure scale height. This assumption is 
based on a simple first order ozone profile determined by 
photochemical equilibrium. Calculations indicate an error 
of less than 5 percent down to 30 km using this approach 
(see Fig. 4-5) . 

4 . 4 Treatment of Multiple Scattering in the Spherical 

Shell Geometry 

4.4.1 The General Problem 

Solutions presented in the previous section (equations 
4-34, 4-35) neglect contributions to the dissociative flux 
due to the scattered or diffuse intensity field at a given 
level. To include multiply scattered radiation in the 
spherical shell geometry, the radiative transfer equation 
(4-5) must be treated with appropriate definition of the 
(^*V) operator and source function J^(r,fi). While the 
problem of multiple scattering in spherical shell atmos- 
pheres has been treated in the study of steller interiors 
(cf. Chandrasekhar, 1960) the problem is somewhat simpler 
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Figure 4-5 

Air mass factors for ozone (Ms oz ) and oxygen absorption and 
Rayleigh scattering (Ms 2 ) . P 




in that case since both the source function and atmosphere 
have spherical symmetry. In the present context however, 
the spherical shell atmosphere is externally illuminated by 
plane-parallel radiation and overall symmetry is lost. 

Thus , while it is relatively straightforward to write down 
the relevant form of the radiative transfer equation (see 
Fymat, 1977) , methods of solution may be tedious. 

An examination of alternative approaches to the prob- 
lem indicates three general methodologies : 

(a) Monte Carlo - (Collins and Wells, 1965; 

Marchuk and Mikhailov, 1967; 

Antyufeyev and Nazaraliyev, 

1973) 

(b) the DART method - (Whitney et al. 1973; Whitney, 

1977; Whitney and Malshow, 

1977) 

(c) approximate semi- 

analytical methods- (Shettle and Green, 1974; 

Sobolev, 1975) . 

Considerations of computational efficiency for radiative 
transfer schemes incorporated within the diurnal photo- 
chemical model severely constrain the applicability of 
approaches (a) and (b) . Therefore, approximate semi- 
analytical techniques have been investigated. They may 
provide viable alternative methods of solution, especially 
when potential errors are assessed relative to other 
photochemical model uncertainties. One such approach is 
the "locally plane-parallel" (LPP) approximation explored 
by Sobolev (1975) . 

4.4.2 Locally Plane-Parallel Approximation 

The conventional plane-parallel (PP) approximation 
assumes a plane-parallel atmosphere uniformly illuminated 
by a plane-parallel source. As an approximate approach to 
treating the general spherical shell problem it may be 
assumed that the atmosphere consists of plane-parallel 
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layers (locally) , illuminated locally, however, accounting 
for the spherical shell nature of the atmosphere (Sobolev, 
1975) . This procedure relaxes one of the assumptions of 
the plane-parallel approximation and should provide results 
with reasonable accuracy. Fundamentally, the methodology 
requires calculations of the primary source function for 
scattering exactly (that is, accounting for Earth curva- 
ture) and approximating multiple scattering using the PP 
approximation . 

By adopting the LPP approach, equation (4-5) becomes 
(wavelength suffix omitted) : 

cos0 yp (t,0) = I(t,0) - J(t,0) (4-36) 


where 0 is the local zenith angle and t(z) is the total 
vertical optical depth to altitude z given by: 


K 

t(z) = E t 1 (z) 
i=l v 


(4-37) 


The source function at altitude t(z) is given by 


t : 


J(t) 


Jo (T) 


(t) 

2 


f +1 

J I (t , y ) dy 


(4-38) 


where y is the cosine of the zenith angle (cos0) 
the single scattering albedo given by: 


U3 0 ( T ) = 


K 


M 11 

a^n /E a n 

k\ • t 
1=1 


0 J o ("r) is 


(4-39) 


This form and subsequent results assume isotropic 
scattering. Anderson and Meier (1979) have demonstrated 
that isotropic scattering provides a good approximation 
to Rayleigh multiple scattering solutions throughout the 
wavelength range .28-. 80 ym. 
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where a D is the Rayleigh (molecular) scattering cross- 
section (Penndorf, 1957), n is the total number density, 
and J q (t) is the primary source function: 

J. (t) = FT [T (z) ,e o 1 (4-40) 


Profiles of optical depth and single scattering albedo vs 
altitude at dissociative wavelengths are illustrated in 
Figures 4-6 and 4-7, respectively. 

The primary source function describes the planetary 
illumination source. In the plane-parallel approximation, 
T is given by equation (4-7) or: 


J 


PP 


(t.) 


0Oq (t) p e — t sec6 o 

4 


(4-41) 


In the LPP approach, the transmission function is defined 
by the spherical air mass factors described in the previous 
section [equation (4-34) 3 and: 

jf P (t) = FT Sp (4—4 2 ) 


4.4.3 Method of Solution 

A variety of techniques are available to solve equa- 
tions (4-36) and (4-38) using the spherical source function 
equation (4-42) in the LLP approximation to obtain the 
intensity field, I, , and ultimately the dissociative flux, 

J ^ 

f“, through equation (4-4) . The approach chosen for this 
application treats the single integral equation of the 
Fredholm type obtainable from equations (4-36) and (4-38) 
by suitable transformation (Busbridge, 1960) : 


J (t) 


Jo (t ) 


. CjOq (t) 
2 



( | T-t ) J ( t) dt 


(4-43) 
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Figure 4-6 

Normal optical depth (t) profile for indicated wavelength. 






where E^_ is the first exponential integral (see Appendix C) 
with general form: 


E (x) 
n 





(4-44) 


and t* is the total normal optical depth [i.e., t* = 
t(z=0)]. When Lambert surface reflection of both the 
direct beam and scattered intensity field is included in 
the formulation, the integral equation for the source 
function becomes : 


J(t) = J 0 (t) + 


CO , 


T 1 l J(t) E i 


( T-t ) dt 


(4-45) 


+ RE 9 ( T * - T ) [2 COS0 o J Q (t* ) + CU o ( T ) / J(t) 

^ o 


E2 ( t * — t ) dt ] 


where R is the surface albedo. The required dissociative 
flux is immediately available noting: 


F d [T (z) ] = 4 


“o [t (z) ] 


J [T (z) ] 


(4-46) 


or alternatively: 

r T * 

F d vx) = FT ( T , 8 0 ) + 2 J J(t)E ( |T-t|)dt 

o 

+ 2 RE 2 (t*-t ) [cosO 0 FT (t* , 0 o ) (4-47) 

+ 2 J J(t)E 9 (T*-t)dt] 

where T[t(z), 0 o ] is the spherical transmission function 
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given by equation (4-34) . The four terms in equation 
(4-47) above represent contributions to the dissociative 
flux due to: (a) attenuation of the direct solar beam, 

(b) the angle integrated diffuse (scattered) intensity 
field, (c) surface reflection of the attenuated direct 
solar beam, and (d) reflection of the surface incident 
angle integrated diffuse intensity field. 

Since the source function J(t), is proportional to the 
single scatter albedo , u ) 0 (t) , [see equations (4-38), and 
(4-40)], the dissociative flux equation (4-47) reduces to 
the absorption only cases [equations (4-6) , and (4-35) ] 
when the albedo is identically zero (i.e., when there is no 

4 - 

scattering ) . 1 

Solution of equation (4-45) is possible using a 
variety of techniques including variational methods (Sze, 

1976) and iterative techniques (Irvine, 1965) . Iterative 
techniques such as successive orders of scattering (SOS) 
are attractive since they easily accommodate the inhomogen- 
eous nature of the atmosphere (as evidenced by the depend- 
ence of single scattering albedo on altitude shown in 

Fig. 4-7) and provide physical insight into the relative 
roles of single vs higher order scattering in determining 
dissociative fluxes. Furthermore, if only a few orders of 
scattering are required to reach a reasonable solution, SOS 
approaches may be computationally competitive with other 
multiple scattering treatments SOS calculations converge 
rapidly to the desired accuracy for small values of single 
scattering albedo since the n*" order scattering contribu- 
tion to the source function is proportional to (Lenoble, 

1977) . Convergence slows as the albedo approaches unity 


t 


With the exception of the surface reflected term which 
was not included in the previous development. 
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(very weak absorption) and computational times become 
prohibitive as total optical depths, t*, become large 
(Nagel et al. 1978) . 

Examination of Fig. 4-7 suggests that the magnitude of 
the single scattering albedo will be small (<0.2) at 
altitudes greater than approximately 15 km for dissociative 

o 

wavelengths <3100 A. Alternatively, when albedos approach 

o 

unity (generally for X > 3400 A) , optical depths are 
generally less than 1.0 (see Fig. 4-6). Based on these 
observations an SOS treatment has been adopted to evaluate 
multiply scattered dissociative fluxes within the wave- 
length region of interest. 


4.4.4 Single Scattering Results 

The first order effects of molecular scattering within 
the spherical shell atmosphere are due to single scatter- 
ing. Since single scattering calculations are relatively 
straight-forward, it is of considerable interest to examine 
the degree of accuracy (compared to multiple scattering) 
obtainable by considering single scattering alone. 

Equation (4-45) is quite amenable to a single scatter- 
ing treatment by replacing the full source function, J, 
appearing in the integral terms with the primary term, J Q 
[equation (4-40)] alone. In the single scattering approxi- 
mation : 


F d (T,T*,e o ,R) = F{T(T,e „)+!,/ f.(t)E 1 ( t-T )dt 

O ■ L 

(4-48) 

r T * 

+ RE 2 (t*-t) [2cos6 o T(t*, 0 o ) + J f D (t)E 2 (T*-t)dt]j 


where 


f D (t) = 0> o (t)T(t,e o ) 


(4-49) 
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and T[t(z), 0 o ] is the transmission function. 

Expression (4-48) is near analytical with the require- 
ment that two integrals be evaluated. An efficient quadra- 
ture technique was developed using piecewise continuous 
functions and exploiting the recursion relation for expon- 
ential ( E n ) functions to evaluate the integral terms in 
equation (4-48) (see Appendix C) . 

Figures 4-8 and 4-9 summarize sample calculations of 

d 

normalized dissociative fluxes (i.e., F /F) at wavelengths 

o o 

of 3475 A and surface albedo 0.2 and 3075 A, R = 0.0, 
respectively, for a solar zenith angle of 60°. Curve (1) 
[dotted] is the result for absorption only, curve (3) is 
the corresponding multiple scattering (MS) result from 
Luther and Gelinas (1976) . Curve (2) is the single 
scattering result consisting of components (a,b,c,d) 
corresponding to the four terms in equation (4-48) above: 

(a) the attenuated direct solar beam, (b) the singly- 
scattered diffuse flux, (c) the surface-reflected direct 
solar beam, and (d) the surface-reflected singly scattered 

O 

diffuse flux. At 3475 A there is virtually complete 
transmission at all altitudes for a purely absorbing atmos- 
phere. Therefore, consideration of scattering has a marked 
effect on the dissociative flux profile. Examining the 
individual contributions to the singly scattered dissocia- 
tive flux (curve 2) in comparison to the MS result (curve 3), 
it is noted that the shape of the MS flux profile is largely 
determined by enhanced extinction of the direct solar beam 
to scattering in the lower atmosphere (curve a) . As noted 
by Lacis and Hansen (19 7 4) , the lower atmosphere back- 
scatters a fraction of incident solar flux (curve b) to the 
upper atmosphere in a manner analogous to an effective 
surface reflection. Contributions due to surface reflection 
(curves c and d) are significant mainly in the lower atmos- 
phere (z < 10 km) . Above 20 km, the SS flux profile 
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NORMALIZED FLUX 

Figure 4-9 

<ji Normalized dissociative flux profile (F^/F) at 3075A, 0 O = 60°, R = 0*0. 



underestimates the MS result by about 16 percent. 

Absorption by the Hartley-Huggins bands of ozone 

o 

dominates the dissociative flux profile for 3075 A at 
altitudes greater than 20 km. Below this level scattering 
contributions modify the profile from the absorption only 
case (curve 1) . The SS results (curve 2) provide a reason- 
able estimate of the altitude dependence but due to 
increased optical depths the MS result (curve 3) is under- 
estimated by as much as 35 percent near the surface. 

Calculations described above are valid in the plane- 
parallel limit (0 O = 60°). The spherical shell nature of 
the atmosphere becomes apparent at zenith angle greater 
than about 80°. Figure 4-10 shows single scattered 
dissociative flux profiles for solar zenith angles of 60° , 
89°, 92° and 96°. At equinox and latitude 30° these 
zenith angles correspond to the local times given in 
Table 4-2 depending on whether dawn or dusk is considered. 
It is noteworthy that in the plane-parallel approximation 
dissociative fluxes at all altitudes are identically zero 
for 0 O >_ 90° (i.e. , t _< 0600 or t _> 1800) . As indicated 
there is no real singularity at dawn or dusk. In the 
approximate 14-minute period between 0 O = 89° and 92°, 
fluxes above 30 km remain virtually unchanged. A feature 
which does appear at 0 O > 90° due to the Earth's curvature 
is the terminator shadow which reaches a height of about 
35 km approximately one-half hour before dawn and after 
dusk. 

In order to assess the relative accuracy of the single 
scattering approach, multi-wavelength calculations of 

O 

dissociative flux profiles in the range 2500-4000 A were 
performed using equation (4-48) for a solar zenith angle of 
60°. Comparisons to multiple scattering results (Luther 
and Gelinas, 1976) for surface albedos of 0.0 and 0.25 are 
illustrated in Figs. 4-11 and 4-12, respectively. Plotted 
are percent errors defined as: 
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ALTITUDE (km) 




TABLE 4-2 


Corresponding Local Times (hr) 


Solar 



Zenith 

Dawn 

Dusk 

Angle 

Hours 

Hours 

60 

0821 

1539 

89 

0605 

1755 

90 

0600 

1800 

92 

0551 

1809 

96 

0532 

1828 
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50 



Figure 4-11 

Percent error profile for single vs multiple scattering 
at indicated wavelength (0 O = 60°, R = 0.0). 
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ALTITUDE (km) 


Figure 4-12 

Percent error profile for single vs multiple scattering 
at indicated wavelength (0 O = 60°, R = 0.25). 







A(%) 


Fmc ( z r 0 o ) - F^ c ( Z , 0 o ) 


x 100% 


(4-50) 



(z , e 0 ) 


It is apparent that accuracies of greater than 25 percent 
can be obtained (for R = 0.25) using a single scattering 
approach at all wavelengths and at altitudes greater than 
approximately 20 km. (Errors are enhanced by increasing 
surface albedo.) This implies that between 70 and 80 
percent of the multiply scattered dissociative flux for 
these cases is due to single scattering alone. Accuracies 

o 

are best for wavelengths (A _< 3300 A) where molecular 
absorption plays. a role since single scattering dominates 
the diffuse flux field. For wavelengths where there is 

O 

little gaseous absorption (A > 3300 A) accuracies increase 
as total optical depth decreases. 

In the lower atmosphere (z < 20 km) where molecular 
scattering dominates (oj 0 ^ 1.0) , errors as high as 40 per- 
cent are encountered at shorter wavelengths. This is due 

to the high attenuation (optical depths > 1.0) of the 
direct solar beam and dominance of the resultant diffuse 
field. Considering other constraints (such as the presence 
of clouds) which may limit more sophisticated approaches in 
the lower atmosphere the accuracy of the single scattering 
treatment may be deemed acceptable. This approach has been 
employed in subsequent calculations of diurnally dependent 
photodissociation rates discussed in §4.5. 

4.4.5 Multiple Scattering Results 

When accuracies greater than those available from 

single scattering calculations are required equation (4-45) 
may be solved by successive orders of scattering (see 
section 4.4.3). This treatment provides additional insight 
into the relative role of higher order scattering in 
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order scatter- 


determining the dissociative flux profile. 

th 

The dissociative flux due to (n + 1) 
ing is evaluated by generalizing equation (4-47) : 


F^ +1 (x) = Fo (T,e o ) +2 / J n (t)E x ( |T-t| )dt 


+ 2R E 0 (t*-t) [cos0 o F o (t*, 0 o ) + 2 


/ 


(4-51) 


J n (t)E 2 (t *- t)dt ] 


where 

Fo(t,0o) = FT (t , 0 o ) 


(4-52) 


and 


J n (x) = F^(t, 0 o ) (4-53) 

For example, the first order (single) scattering solution 
[equation (4-48)] corresponds to n = 0 . 

Trial calculations based on equation (4-51) are illus- 

o 

trated in Fig. 4-13 for a wavelength of 3475 A, solar 
zenith angle of 60° and surface albedo of 0.25. At this 
wavelength the single scattering albedo profile is near 
unity throughout the atmosphere (Fig. 4-7) and the total 
optical depth is about 0.6 (Fig. 4-6). Based on these 
optical parameters, calculations in this wavelength region 
represent a worse case situation for multiple scattering, 
i.e., one where the contributions to the dissociative flux 
due to multiple scattering are maximum, especially in the 
lower atmosphere. It is instructive, therefore, to ascer- 
tain the order of scattering required to achieve the 
multiple scattering result for this constraining case and 
the improved accuracy (with respect to multiple scattering) 
incurred with each successive order of scattering. Conver- 
gence to the multiple scattering solution (within 0.5%) 
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Figure 4-13 

Multiply scattered dissociative fluxes by successive 
orders treatment where n is the order of scattering 
(A = 3475A, 0 O = 60°, R = 0.25). 
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requires 9 iterations, however, even a second or third 
iteration drastically improves the result. For the case 
illustrated as much as 95 percent of the dissociative flux 
in the upper atmosphere may be attributed to three orders 
of scattering. 

This may be seen explicitly in Fig. 4-14 which illus- 
trates the percent error defined as : 


A (%) 




x 100% 


(4-54) 


as a function of iteration number for a wavelength of 

O 

3275 A at altitudes of 10, 25 and 40 km. To achieve 
comparable levels of accuracy approximately one more 
iteration is required in the lower atmosphere (10 km) than 
in the upper atmosphere. In this example, the single 
scattering albedo in the lower atmosphere is approximately 
unity while the mean albedo in the upper atmosphere 
(z > 20 km) is about 0.5. The total optical depth is 
closer to unity. As expected based on the previous discus- 
sion of the SOS treatment for multiple scattering, conver- 
gence to a desired accuracy requires consideration of fewer 
orders of scattering for albedos less than unity (i.e., 
above 20 km) . 

With the methodology for treating multiple scattering 
formulated in this manner (as a logical extension of the 
single scattering approach) and the required algorithms in 
place, comprehensive sensitivity analyses of wavelength 
dependent dissociative fluxes and photodissociation rates 
to higher orders of scattering may be undertaken. However, 
computational costs for each wavelength increase linearally 
with the order of scattering considered. Due to the finite 
resources available within the present study, the decision 
was made to concentrate, therefore, on constraining cases 
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ORDER OF SCATTERING 


Figure 4-14 

Percent error vs multiple scattering result for indicate^ 
order of scattering at altitudes of 10, 25, 40km (X=3475A, 
0 o = 60° , R= 0 . 25) . 
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such as those described above. For similar reasons, 
results presented in subsequent sections are based on 
terminating the multiple scattering code after one itera- 
tion. As previously discussed, this insures accuracies in 
dissociative fluxes on the order of 20-30 percent through- 
out most of the atmosphere above the tropopause. 

4 . 5 Application to Diurnally Dependent Photodissociation 

Rates 

As a prerequisite to calculation of diurnal variations 
of radical species in the spherical shell geometry, techni- 
ques described above were applied to evaluation of photo- 
dissociation rates. Single-scattered dissociative fluxes 
[equation (4-48)] were computed as a function of wave- 
length, altitude, and solar zenith angle (time of day) and 
corresponding dissociation rates were evaluated using 
equation (4-3) and appropriate cross-section data. Table 
4-3 lists the set of dissociative reactions treated in the 
analysis. For example. Fig. 4-15 illustrates the vertical 
profile of the ratio of single scattered (SS) dissociation 
rate to that for pure absorption for reaction 11 in Table 
4-3 evaluated at a solar zenith angle of 45° and surface 
albedo of 0.2. For comparison, multiple scattering (MS) 
results are presented obtained using a variational approach 
to solve equation (4-45) (Sze, 1976). Multiple scattering 
results for a slightly higher surface albedo (R = 0.25) are 
also included (Luther et al. 1978). Results indicate that 
single scattering accounts for 75-85 percent of the total 
(i.e., MS) photodissociation rate in the upper atmosphere 
(z > 20 km) . 

The difference between plane-parallel and spherical 
shell geometries is highlighted by examining variation of 
computed photodissociation rates at times near the local 
terminator. Figure 4-16 shows the rates for reactions 11 
and 20 in Table 4-3 at altitudes of 25 and 40 km near dusk. 


86 



TABLE 4-3 


Photolysis Processes 
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Figure 4-15 

Comparison of photodissociation rate for process N02^ NO+O evaluated by successive 
order treatment truncated after one iteration (i.e., single scattering) to multiple 
scattering result (0 O = 60°, R = 0.2). 









As indicated in the inset, sunset at the surface occurs at 
1800 for the case considered (equinox, 30° latitude) . In 
the plane-parallel treatment, photodissociation rates at 
all altitudes become identically zero at this time. As 
indicated by the terminator position times for 25 and 40 km 
(1823.4 and 1829.6, respectively), the first order effect 
of sphericity at elevations above the surface is to extend 
the length of the day . 1 This amounts to about 3/4 and 1 
hour, respectively, at 25 and 40 km (due to dawn/dusk 
symmetry) . 

4 . 6 Diurnal Calculations 

Time dependent altitude profiles of free radical 
species (e.g. , CIO, C 10 N 02 * HC1, H 2 O 2 , HO 2 , OH, NO, NO 2 ) 
were obtained by performing diurnal calculations using 
numerical procedures described in Appendix D. Photodissoc- 
iation rates used in the diurnal simulations included: 

(a) molecular scattering effects, (b) surface reflection 
(a surface albedo of 0.2 was used throughout the wavelength 
region) and (c) consideration of the spherical shell geo- 
metry of the Earth's atmosphere. Treatment of these 
factors has been discussed in previous sections. Single 
scattered dissociative fluxes were employed as a first 
order treated. As noted above (§4.5) this approximation 
accounts for the bulk of molecular scattering within the 
wavelength interval. In conducting the diurnal calcula- 
tions, long-lived species (such as N 2 O, O 3 , CH 4 , etc.) were 
held fixed and vertical transport was neglected. 

Figures 4-17a and b illustrate time dependent results 
for CIONO 2 , OH, NO 2 , CIO and HO 2 . An altitude of 40 km was 
chosen for demonstration purposes. The solid species 


+ 


These times neglect refraction. 
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Figure 4-17b 

Comparison of diurnal variation of NO 2 / CIO, and HO 2 using 
absorption only-plane parallel treatment (dotted) and 
scattering-spherical treatment (full) at z = 40km. 
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abundance curves refer to the spherical shell (Sp) , single 

scattering result. For comparison, plane-parallel (PP) , 

pure absorption results (dotted curves) are also provided. 

The time scale refers to local time at equinox and 30° 

north latitude. Times of dawn (0600) and dusk (1800) at 

+ 

the surface (z-o km) are indicated by light vertical 
lines . 

4.6.1 Role of Sphericity 

Comparing the two sets of results, the gross effect of 
sphericity on the diurnal variation of free radical abund- 
ances is to accelerate the onset of photochemical produc- 
tion (of OH, CIO and HO 2 ) and loss (of CIONO 2 and NO 2 ) at 
dawn and delay termination of these processes at dusk by 
approximately one-half hour. Thus, at 40 km sphericity 
effectively increases the length of the day by about an 
hour (neglecting refraction) . Since many photolysis 
processes have time constants of seconds to minutes, 
specification of precise local dawn and dusk (i.e., at 
altitude) is of considerable interest for determining the 
transition times between daytime and nighttime chemistries. 

As a consequence of these "phase" lags, the plane- 
parallel approach may seriously miscalculate abundances of 
photochemically active species in the dawn/dusk transition 
period. A case in point is the dawn OH concentration 
(Fig. 4-17a) . In the plane-parallel treatment this is a 
nighttime value, however, when Earth curvature is consider- 
ed the atmosphere at this time is already about one-half 
into daylight and the corresponding OH concentration is two 
orders of magnitude higher. The result is not as striking 
for CIO (Fig. 4-17b) whose major source (photolysis of 


t 


Or at all altitudes in the plane-parallel approximation. 
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CIONO 2 ) has a longer lifetime. 

Several gross features of the diurnal variation of 
those species illustrated in the figures are worth noting. 
These include: (a) the marked difference between daytime 

and nighttime species abundances alluded to above (some- 
times varying by orders of magnitude) , (b) the relatively 

slow variation of concentrations with time during mid-day 
and nighttime hours, (c) the contrasting rapid variation 
during the post-dawn/dusk transition periods and (d) the 
absence of symmetry about local noon (1200 hr) . 

These properties may be understood by considering the 
dependence of modeled abundances on the presence of key 
photolysis processes during the day, their absence at 
night, and the rapid variation of the primary source of 
dissociative flux (the Sun) during dawn/dusk transitionary 
periods. The relative constancy of concentrations near 
midday validates the uses of steady state simulation 
approaches for many applications. One aspect which may not 
be immediately intuitive is the daily asymmetry with 
respect to noontime values. Although photodisspciation 
rates (which to first order are functions of solar zenith 
angle only) are symmetric about noon, photolysis processes 
must compete with other production and loss mechanisms 
which do not necessarily "follow the Sun". Post-dawn and 
post-dusk time dependence is characteristic of the locally 
dominant chemistry and thus, distinctive time constants 
contribute to the post-dusk apparent log linear behavior. 
The attributes discussed above with respect to Figures 
4-17a,b are evidently dependent on the particular species 
under discussion, however, additional calculations indicate 
that they also depend on altitude. 

4.6.2 Role of Molecular Scattering 

The near noontime values shown in Figs. 4-17a and b 


94 



provide a calibration point on the relative role of 
scattering. For these times of day, the present calcula- 
tions reduce to the plane-parallel case and the only 
difference between solid and dotted curves is the treatment 
of first order scattering in the former case. As illus- 
trated, the difference between the two calculations is not 
appreciable. However, the magnitude and sign of the con- 
centration difference near noon is consistent with other 
calculations of the effect of multiple scattering (Luther 
et al. 1978) . At this altitude (40 km) consideration of 
scattering appears to be of minor consequence, however, 
additional calculations indicate that the effect is 
appreciable at lower altitudes. This behavior is under- 
standable upon examination of individual contributions to 
the relevant dissociative fluxes (see Figs. 4-8, 4-9; 
§4.4.4). They suggest that the attenuated direct solar 
beam dominates in the upper atmosphere while diffuse 
contributions become progressively more important toward 
the troposphere (i.e., within 2-3 pressure scale heights of 
the surface) . 

4 . 7 Implications for Remote Sensing 

Previous sections have described a methodology for 
including the effects of Earth curvature and molecular 
scattering in the simulation of diurnally dependent free 
radical species concentrations. Ultimately, the utility 
of extensions to existing theory may be evaluated based on 
their application to either the interpretation of new 
measurements or existing data inconsistent with previous 
theory. In this section implications of the work described 
above are discussed with reference to emerging measurement 
techniques based on remote sensing. 

Historically, most relevant species profile data have 
been collected during daylight hours. This choice is 
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based in part on logistical constraints and in part to take 
advantage of the relative constancy of species concentra- 
tions during mid-day. Such data sets obviate the need for 
other than steady state modeling for most practical 
purposes . 

Requirements for global monitoring of stratospheric 
trace species have motivated recent developments in remote 
sensing technology including use of the solar occultation 
approach (Russell and Drayson, 1972; Russell et al . , 1977; 
Chu and McCormick, 1979; Menzies, 1979). An inherent 
characteristic of the experimental configuration in this 
technique is a sensor field of view with line of sight 
integrating across the terminator. As illustrated in the 
previous section, the number densities of photochemically 
active species may vary by orders of magnitude during dawn 
and dusk. It is evident from these results that since 
there is such a radical time rate of change, photochemical 
equilibrium (i.e., steady state) models may introduce 
uncertainties in the interpretation of data obtained by 
such techniques. 

In principal, most inversion methods for obtaining 
vertical species profiles from occultation data such as 
analytical treatments (Roble and Hays, 1972) and the onion 
peel approach (McKee et al., 1969) rely on the assumption 
of spherical symmetry and stratification in the distribu- 
tion of the relevant absorbers. Since lines of sight in 
the occultation mode essentially look through varying local 
times and altitudes, this requires that observed species 
are time invariants (i.e., depend on altitude only). The 
gross features of the diurnal calculations described in 
§5.6 suggest that this is not true for many radical 
species. Qualitatively, it is expected that inversion 
results based on spherical symmetry considerations will be 
least affected for individual species with suppressed 
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diurnal variations (such as CIO as contrasted with OH or 
NO) . Furthermore, the character of the dawn/dusk asymme- 
tries noted previously suggest that, due to the generally 
more rapid variation of photochemically active species in 
the post-dawn transition period, reduced accuracies may be 
expected for dawn occultation as compared to dusk occupa- 
tions. These assessments are generally confirmed by recent 
simulations of inversions (Boughner et al . , 1980). 

Based on these considerations it may be useful to 
investigate approaches to incorporating a priori informa- 
tion on diurnal variations within the inversion algorithm. 
In this context the occultation techniques can be outlined 
as follows. The radiation of wavelength X received at the 
sensor for a line of sight s with tangent height z Q (see 
Fig. 4-1) will be given by the solution to (4-36) : 

s (z Q ) 

I x [s(z„)] = I x (~)T x [s(z„) ] + J J x (s , )dT x (s') (4-55) 

where 1^ (°°) is the unattenuated solar radiance, T-^ is the 
transmission function along the line of sight equation 
(4-8) , and J is the local source function [determinable 
from equation (4-45)]. If scattering can be ignored, this 
reduces to: 

I[s(z 0 )J = I(co) T [s (z 0 ) ] (4-56) 

where [in analogy to equation (4-9) ] : 

K . r s( ? o) 

T = exp - E a 1 J n 1 (s I )ds l (4-57) 

i=l 

Ignoring refraction, the line integral can be written in 
terms of altitude noting from the geometry: 
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ds (z D ,z ' ) 


(4-58) 


= ( z 1 + R) [ z 1 +R) 2 + ( z 0 + R) 2 ] ^dz 1 
= G(z o; z ' ) dz ' 

yielding 

[ I K . , oo , 

2E a 1 J n 1 (z , )G(z 0 z')dz' (4-59) 

This expression assumes that the relevant absorber is a 
function of altitude only (i.e, spherically symmetric) and 
is amenable to inversion by one of the techniques described 
above . 

Diurnal simulation provide an understanding of the 
true variation of photochemically active absorbers along 
the line of sight. Figs. 4-17a,b, for example, essentially 
indicate : 

n(s') = n[z',0(t)] (4-60) 

where 6 is the local solar zenith angle given by: 

8 = sin ^"[(R + z o ) / (R + z')] (4-61a) 

on the solar side of the path and 

0 = tt - sin "*"[(R + z 0 )/(R + z 1 )] (4-61b) 

on the anti-solar side where 0 = tt/2 denotes the terminator. 

Noting diurnal variations equation (4-59) may be general- 
ized as: 


t (zj 


K 

£ o' 

i=l 


I 


[n (z ' , 0 ) + 


n 1 ( z 


tt- 0 ) ] G (z 0 , z ' ) dz ' 


(4-62) 
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A form analogous to equation (4-59) is obtainable by intro- 
ducing a shape function describing the spherical asymmetry 
obtainable from diurnal simulations (Boughner et al., 1980): 


D (z 0 , z 1 ) = 


r n^(z , f 6) + n~*~ (z 1 ,tt- 6 ) 
2n 1 ( z ’ , tt/2 ) 


yielding : 


t (z 0 ) 


K . r 

2E o 1 J n 1 ( z 1 , n/2 ) K 1 (z 0 ,z 1 ) dz 1 
i=l z 0 


(4-63) 


(4-64) 


where the kernal is : 

K 1 (z 0/ z f ) = D 1 (z 0 ,z , )G(z 0 ,z l ) (4-6,5) 

For time invariant species D 1 = 1 and equation (4-64) 
reduces to the spherically symmetric case. 

Equation (4-64) is a Fredholm equation of the first 
kind defined generally as: 

r b 

h (y ) = J K(x,y)f(x)dx (4-66) 

a 

where h (y) represents the data and f (x) the desired 
function. In operator form equation (4-66) may be written 
as : 

Kf = g (4-67) 

In this form a number of powerful generalized inverse tech- 
niques such as Phillips (1962) and Twomey (1950) are appli- 
cable which do not assume spherical symmetry. Furthermore, 
questions regarding uncertainties introduced in simulating 
the required kernal function equation (4-63) and other 
typical sources of noise such as those introduced in 
linearization and quadrature of equation (4-64) may be 
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addressed by examining optimal smoothing parameters, 
selecting suitable regularization functions, and employing 
appropriate inversion constraints (such as non-negativity) . 
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5. 


SUMMARY AND RECOMMENDATIONS 


This section summarizes some of the major findings of 
the tasks described in the previous sections. The reader 
is referred to the individual sections for more detailed 
description of results. 

5 . 1 Model Sensitivity Studies 

As discussed in Section 2 the results of 1-D model 
simulated ozone perturbations to high altitude aircraft 
operations are subject to major uncertainties particularly 
as regards rate constant data related to OH chemistry. 
Sensitivity studies indicate that calculated OH concentra- 
tion profiles are most sensitive to values adopted for the 
reactions listed in Table 2-1. Perturbation calculations 
were performed for three model scenarios corresponding to 
high (model A) and low (models B and C) values of strato- 
spheric OH concentration. Model C investigates the poten- 
tial role of peroxynitric acid (HO 2 NO 2 ) as a sink for OH 
and provides indicators of OH concentration (e.g., CIO , 
HNO 3 /NO 2 r HF/HC1) in best agreement with measurements. 

Model A tends to predict an increase in column ozone for 
both low (15 km) and high (20 km) NO x injections. The low 
OH models (B and C) , however, predict a small (<1%) increase 
for the low altitude injection, but a fairly significant 
reduction for the high altitude case (-5.5%, model C) . More 
reliable kinetic data is needed to differentiate between 
the models discussed. 

Additionally, it is noted that injection of NO x may 
perturb both local stratospheric and surface temperatures. 
Thus, climatic changes should be considered in addition to 
column ozone perturbations. 
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5 . 2 Two-Dimensional Zonal-Mean Modeling 


The application of Eulerian averaging processes to 
zonal-mean models introduces a number of undesirable 
features associated with the eddy flux contributions to 
transport. The feasibility of some kind of Lagrangian 
approach as an alternative averaging methodology is attrac- 
tive from the standpoint that eddy terms are avoided. 
However, difficulties arise in the implementation of the 
formulation and in interpretation of mixing ratio results 
in comparison with observations. Several areas are identi- 
fied for additional research before the Generalized Lagran- 
gian Mean (GLM) tracer equation may be applied to zonal 
modeling . 

In the short term, numerical tracer experiments compar 
ing Eulerian and GLM approaches may be performed deriving 
the required Lagrangian velocity components for given three 
dimensional motion fields. One potential avenue for 
research would concentrate on cases with small disturbances 
using a more precise treatment of the diabatic heating rate 
to generate Lagrangian velocity fields. 

Longer term research may be profitably directed to 
examining the interpretation of the basic continuity equa- 
tion and in reformulating the momentum equation by defining 
more generalized differential operators. 

5 . 3 Remote Sensing 

The role of multiple scattering and Earth sphericity 
on the calculation of diurnally dependent photodissociation 
rates and trace species concentrations was examined in 
Section 4 with particular emphasis on dawn and dusk. 

Calculations of slant path air mass factors indicated 
that sphericity should be considered in evaluating the 
primary source function for scattering for solar zenith 
angles greater than about 80°. Evaluation of profiles of 
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atmospheric optical depth and single scatter albedo at dis- 
sociative wavelengths suggested that a successive order of 
scattering (SOS) approach may be viable at all altitudes. 

A locally plane-parallel (LPP) treatment of the radiative 
transfer equation was adopted to treat the spherical geom- 
etry. Comparison of single scattered dissociative flux 
profiles with multiple scattering results indicated that 
between 70 and 80 percent of multiply scattered flux above 
20 km can be accounted for by single scattering alone. In 
the lower atmosphere where molecular scattering dominates 
errors as high as 40 percent are encountered, however, 
accuracies comparable to those in the upper atmosphere can 
be achieved by considering an additional order of scatter- 
ing. Convergence to the multiple scattering result (within 
0.5%) requires about 9 orders of scattering, however, 

5 percent accuracy is available in the upper (lower) atmos- 
phere with only 3 (4) orders. Diurnal calculations of 

trace species indicate that the gross effect of sphericity 
on diurnal variation is to accelerate the onset of photoly- 
tic processes at "plane-parallel" dawn and delay their 
termination at "plane-parallel" dusk. 

Based on the results of diurnal trace species varia- 
tion across the terminator region, it is recommended that a 
potential approach to the inversion problem in the occulta- 
tion geometry may be consideration of generalized inverse 
techniques not relying on assumptions of spherical symmetry. 
In this manner the shape of the kernel function may be 
estimated to first order using a priori calculations of 
diurnal trace species concentrations. 
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APPENDIX A 


Earth-Sun Geometry 


As in the plane-parallel case, the transmission path 
for solar radiation incident at a specified altitude, z(km), 
will be a function of the local solar zenith angle, 0 O . 

The solar zenith angle may be defined from the scalar 
product of the local surface normal vector, n, and the 
Sun's direction vector, s. In a celestial coordinate 
system with origin fixed at the Earth's center, the unit 
normal vector will be a function of latitude, A, and local 
hour angle, HA. The unit solar direction vector, s, will 
depend on the Sun's declination angle, 6. The required 
relationship between solar zenith angle, latitude, solar 
declination angle and hour angle is given by: 

cos0 o - sin6 sinA + cos6 cosA cosHA (A-l) 

Solar declination angle, 6, may be approximated as: 
sin0 = 0.3978 sin [ . 98 63 (d-80 ) ] (A-2) 


where d is the day of the year {1 Jan. = 1, 20 Mar. = 80, 
etc.) . Note that 6 is zero at the equinoxes, and has 
maxima 1 6 f = 23° 27' at the solstices. 

The hour angle, HA, is defined as identically zero at 
local noon, 1200h, (Sun's direction vector along observer's 
meridian) when the Sun's zenith angle is smallest: 

cos6 o (1200) = sin6 sinA + cos6 cosA (A-3) 


Conversely, sunrise and sunset times may be defined as 
corresponding to solar zenith angles of 90°. The length of 
the day in hours is then given by: 


LOD 


2 -1 

cos [-tan6 tanA ] 

ID 


(A-4) 
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Thus hour angle may be written as : 


HA (deg) = 15(12.00-t) + e (A-5) 

where 15 is the Earth's angular rotation rate (°h and e 
is a correction for the equation of time. Equations (A-l) , 
(A-2) and (A-5) define a solar zenith angle for a given 
date, time, and latitude. 

Additionally, the local normal, n, and the Sun's 

A 

direction vector, s, define a plane of incidence pictured 
in Fig. (A-l) . This plane contains the subsolar point (SSP) 
where the solar zenith angle is zero and the required local 
normal vector, n. For solar zenith angles which are 
sufficiently small (say 0 O < 80°), a plane-parallel treat- 
ment will suffice. For situations near sunrise and sunset, 
however, the Earth's sphericity is a critical consideration. 
Figure (A-l) illustrates a number of relevant cases. Local 
normal n denotes a situation with 0 O approaching 90°. In 
this case the usual plane parallel treatment (sec0 o ) will 
overestimate the relevant optical path. As illustrated by 
local normal r\2 / the solar zenith angle 0„ may be >90° at 
elevations above the surface. For these situations the 
upper atmosphere may be illuminated while the lower atmos- 
phere is not. The shadow height will be given by: 


Z g = R(csc0 o - 1) 
= R(secU - 1) 


(A-6) 


where U is the solar depression angle defined as: 

0o (>90°) = 90° + U (A-7) 

and R is the Earth's radius ( 6371x10 3 m) . Thus for U - 8° 
sunrise occurs at an altitude of 62.6xl0^m. This corres- 
ponds to about 0523 local time at equinox at 30 °N. [These 
calculations and equation (A-6) ignore refractive effects 

which decrease Z . ] 

s 
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APPENDIX B 


Procedure for Refracted Air Mass Calculation 


Calculations follow a procedure based on Selby and 
McClatchey (1972). Essentially, equation (4-22) is modi- 
fied by a factor dependent on the phase refractive index of 
air, m(z) with the resultant air mass factor given by: 



(z,0 o ) 


1 

n 1 (z)H z 



(z') 


( R+z ) m ( z ) 2 

( R+ z ' ) m ( z ' ) 


■ 2 q 
sm 0 



dz ' 


(B-l) 


The phase refractive index of air m(z) is given by Edlen's 
(1966) expression: 


[m ( z ) - 1.0] x 10° = (77.6 + 


0459 P (z) 


T(z) 


p h 2 o ^ 43 * 49 


0.347 

A 2 


(B-2) 


where A 


P 

T 

p h 2 o 


wavelength in microns (vim=10^nm) 
pressure (atm = 101.325 kPa) 
temperature (K) 

partial pressure of water vapor (atm = 

101.325 kPa) 


The geometry of the refraction calculation is illus- 
trated in Fig. (B-l) . Following Selby and McClatchey 
(1972) , the path integral is performed using the local 
solar zenith angle as the integration variable. For 
initial solar zenith angle + ^ incident on level Zj + ^, 

geometry yields: 
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Figure B-l 

Geometry of refraction calculation. 
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(B-3) 


RH ' Z 1 + 1 
sin0 j 


R+z . 

i_ 

sina . 
3 


From Snell's law: 


sma. 


mj +! . 

IT sine j+i 


(B-4) 


The incidence angle on the next level is obtainable from 
equations (B-3) and (B-4) : 


sinQ j 


( R+z j-n) frj+i 


(R+Zj) 


m. 


sin6 j+l 


(B-5) 


or 


m sin0 (R+z) = const 


(B-6) 


The increment of path length in the layer from z^ to 
is given by: 


As 


j + 1 


(R+z 



sing . 

1 

sin8 . 
3 


(B-7 ) 


where 


is given by: 


0 . 
3 


a . 
3 


(B-8) 


Substitution of equations (B-4) , (B-5) and (B-8) into (B-7) 

yields : 


As j+ i 


< R+z i + i> 
sin0 • 

3 


. -1 


sm \ sm 


f ( R +z j +1 ) m j+i ^ 

(R+z.l ^ Sln6 j+lj 


. -1 
- sm 


m 




m. 


sin0 . , 
3 + 1 


(B-9) 
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APPENDIX C 


Integral Quadrature Technique 

Evaluation of expression (4-48) for single scattering 
and the general successive order of scattering relation 
(4-31) require calculation of integrals of the form: 


f T 

= J E ( |t — t | ) f(t) dt 

II 1 1 


(C-1) 


th 

where the n exponential integral is given by: 


n 


/ \ f - x/y n-2 

(x) = J e / H y 


dy 


(C-2 ) 


The first three exponential integrals are illustrated in 
Figure C-1. These integrals may be approximated as: 


N r 

I (t , t*) - Z I J X 

n 'v . J 

1 = 1 t . 

J 


E n ( |t-T | ) f (t) dt 


(C-3 ) 


where t^ = 0 , t 2 = • . . , t 


N+l 


= T* 


Defining : 


F j (T) 


t 

= J ^ +1 E n (lt-i|) f(t) dt 


t . 
J 


( C— 4 ) 


yields 


N 


n 


I (t ,t*) - Z F . (t) 


n 


* j=i 3 


(C-5) 


The appropriate quadrature formula is of the form: 


F^(t) = a^(T) fj + a" + 1 (T)f j+1 


(C-6) 
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Figure C-l 

Exponential integral functions E-^ , E2 / E3 (x) . 


Ill 





where : 


(C-7) 


fj - f(tj) ? f j +1 

If piecewise continuous (local linear interpolation) 
tions are used: 

(t-t . ) 


f (t> = fjU - + f. 


(t-tj) 


j+1 At.. 


where 


At . = t . , - t . 

3 1 + 1 3 


Substituting into equation (04) : 

t_; 
t 


t . , 

'j(T) = fj J 3 E n ( It-T I) [1 - 


(t-t . ) 
— ] dt 


At. 


r (t-t.) 

+ f i+l /. E n ( |t-x|)-^_ dt 


t . 
3 


From equation (C-6) : 


r t . , , (t-t . ) 

i = I E n ( 11 At 1-1 ^ 


t . 
3 


n 

a j+l 


f t • ^ 


Note that : 


• + *j +1 = / tj+1 E n ( |t-T,)dt 


func- 
(08 ) 

(09) 

(C-10) 

(011) 

( 012 ) 

(013) 
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Expressions (C-ll) and (012) may be evaluated analytically 
by using the recurrence relation: 


E n (x) dx = -d E n+1 (x) 


(C-14 ) 


Taking note of absolute values the constants a 1 ?, a 1 ?,, 

3 3+1 

obtained are : 


a j (T) = E n+l (t j- T) 


E . 0 (t • “T ) - E 


At j \ n+2 v j 


n+2 


(t j+ i- T )} 


(C-15) 


n 


a- (t) = -E„. 1 (t. ,,-t) + 


j + 1 


n+1 j + 1 


_J_ / 

At. \ 


E n + 2 (t j' T) 


E n+2 (t j+l T) } 


(C-16) 


for t < t . and 
- D 


n / v 
a j ( t) = 


-E (t 
n+1 


“V - aft ! 1 


E n+2 (T E n+2 (T ‘ t j+l ) 


(017) 


n / v 
a j+i T 


“ E n+l (T - t j+l> + At? ^ E n+2 (T_t j ) 


E n+2 (t 


t j+l ) ) 


(C-18) 


for t > t . , , 
- D + l 
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APPENDIX D 

Numerical Procedures for Diurnal Calculation 


The distribution of an atmospheric trace gas is 
governed by the continuity equation, 

~(i) + ||(i) = P(i) - L ( i) n (i ) (D-l) 

and the flux equation 


<f> (i) 



( 


n (i) x 
N ; 


( D — 2 ) 


th 

where n(i) is the number density of the i species, (J) (i) 
the vertical flux, P(i) and L(i)n(i) the chemical produc- 
tion and loss terms, K the eddy diffusion coefficient and N 
is the total air number density. In general, the chemical 
production and loss terms depend on other species. Equa- 
tions (D-l) and (D-2 ) represent a system of time dependent, 
coupled and nonlinear differential equations with a wide 

_ O f. 

range of time scales (10 -10 sec) appropriate for various 

species of interest. 

For many applications, such as the interconversion of 
active chemical species, the transport term, 3<J)/9z, is 
orders of magnitude smaller than the chemical terms and may 
therefore be neglected. Thus, equation (D-l) reduces to 


dn (i) 
dt 


P (i) - L (i) n (i) 


(D-l') 


subject to periodic boundary conditions [i.e., n(i, t=0) = 
n(i, t=T) , where T = 1 day]. 

Equation (D-l) is replaced by a finite difference 
equation (Wofsy, 1977), 

n £ +1 (i) " n d (i) " t P l + l (i) - L S,+l (i)n )l + l (i)] At (D " 3) 
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where l is the time index. 

Equation (D-3) represents a fully-implicit scheme and 
must be solved by an interactive algorithm at each time 
step. It can be easily shown that the finite difference 
scheme equation (D-3) has the following important proper- 
ties : 

1 ) satisfying mass conservation at every time step; 

2) stable for both P and L at every time step. 

Equation (d- 3) is used to solve for the concentrations 
Of 0( 3 P), 0( 1 D), 0 3 , NO, N0 2 , HNO 3 , N0 3 , N 2 0s, Cl, CIO, 

HC1 , C10N0 2 , H, OH, HO 2 r H 2 0 2 . For a species with chemical 
time constant shorter than a day, only a few days are 
required to achieve periodic boundary conditions. However, 
for species with chemical times significantly longer than 
one day, over a hundred model days are usually required to 
achieve periodic boundary conditions. We may remedy this 
situation by setting the initial conditions 

r t r T 

i o P ( i , t ) exp [ - / L (i,s ' ) ds ' ]ds 

n ( i , t =0 ) = - — ( D-4 ) 


1 - exp[ - f T L (i , s ' ) ds '] 
J O 


where T = 1 day. Clearly, for long-lived species the loss 
frequency L(i,t) satisfying 


J\a, 


l , t) dt << 1 


(D-5 ) 


Equation (D-4) reduces to 

/ T P(i,t)dt 

n(i, t= 0 ) = 


[ L ( i , t ) dt 
J O 

In fact, if equation (D-5) is satisfied, then we have 
n (i , t) = n (i , t=0) for 0 < t < T 


(D-4') 
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This procedure of adjusting initial conditions at the 
beginning of each model day greatly reduces the number of 
model days required for relaxation. 
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APPENDIX E 


Initialization of Langrangian Trajectory 

[The content of this Appendix is based on work from 
Tung, private communication.] 

o 

Let £' (t,x)be the particle displacement trajectory 

o 

of the particles located at coordinates x at time t Q . It 
follows that the position of the particle at subsequent 
time t is given by 

S(t,x) = x + |* (t,x) (E-l) 

We define the coordinate x 0 , and £(t,x 0 ) by 


x Q — S ( t , x) ( E-2) 


| (t ,x 0 ) = S ( t , x ) - x 0 


(E-3) 


Note that equation (E-3) implies 


£ (t ,x 0 ) = 0 


(E-4) 


If in addition £(t,x 0 ) satisfies equation (3.3.2b) , £ can 
then be used as a disturbance induced displacement vector 
in Lagrangian averaging. 

The Eulerian velocity v(t,S) is given by 

^•S = v (t,S) (E-5) 

From equation (E-3) it follows that 
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v(t,S) = ^x 0 +^|(t,x 0 ) (E-6) 

S 3X 0 

= v(t,x 0 ) + -^£(t,x 0 ) + V x 5(t,x 0 ) 

= v(t,Xo) + (-^ + v(t,x 0 )*V x )|(t f x D ) 

If we apply the Eulerian mean operator to equation (E-6) we 
obtain 

v (t ,S) = v(t,x Q ) , ( E-7 ) 

since x Q = x Q -»■ v(t,x D ) = v(t,x 0 ) and £(t,x c ) = 0. 

However, from the definition of the Lagrangian mean, we 
have 

v (t,S) = v(t,x 0 + £(t,x 0 ) = v L (E-8) 

Combining equations (E-7) and (E-8), we obtain 

v L = v(t,x c ) t (E-9) 

Substituting equation (E-9) into equation (E-6) , we have 

(3^ + Y L * v x o ) | (t ^o ) = “ Y L ^ t,x °^ (e- 10 ) 

^The result here is similar to that presented in Andrew and 
McIntyre (1978) pp. 615-616. 
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which is identical to equation (3.3.2b). Thus, £(t,x 0 ) 
can be used for Lagrangian averaging and is directly 
related to trajectory motion of the particles. 

Note that in the above definition, the GLM operations 
depend on the choice of time t 0 at which the particles are 

O 

labelled by x. In terms of £(t,x 0 ), this is characterized 
by the fact that £(t 0 *x 0 ) =0. A different initialization, 

. . ~ ~ o f 

i.e., choice of t$ and x* can lead to different averaged 
values for some quantities (see McIntyre, 1979). 
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APPENDIX F 


High Altitude Aircraft Emissions 


As a consequence of projected high altitude aircraft 
operations, appreciable quantities of combustion products 
of fossil-origin fuels may be emitted into the ambient 
atmosphere. Typical aviation fuels consist primarily of 
a mixture of higher hydrocarbon petroleum distillates (e.g.. 
Jet A - "kerosene', C9 H 20“ C 16 H 34 ; Jet B - "naptha", 
c 6 h 14“ c 7 h 16^ and inherent residual impurities such as 
sulfur and heavy metals . Exhaust product constituents may 
be categorized as: (1) products of complete combustion 

including water vapor (H 2 O) and carbon dioxide (CO 2 ) , 

(2) products of incomplete combustion including carbon 
monoxide (CO) and hydrocarbons (HC) , and (3) products due 
to residual impurities such as sulfur (S) and its oxides 
(SO 2 ) #■ and oxides of nitrogen (N0 X ) . This discussion 
focuses on the evaluation of aircraft related yearly source 
strengths (MT yr 1 ) for H 2 0, CO, S and N0 X . Total emissions 
are dependent both on intrinsic emission characteristics of 
individual aircraft/engine configurations as given by the 
emission index and fleet related modes of operation includ- 
ing total number of aircraft and cruising time per day. 

Emissions of individual pollutants and pollutant 
groups are quantified in terms of relavant emission indices, 
El, which are the number of grams of pollutant produced for 
each kilogram of fuel consumed. In general, El are depend- 
ent on engine type, mode of operation (cruise, idle, etc.), 
cruise power setting, flight altitude, and cruising speed. 
However, El have been recommended for typical cruise 
conditions for performing atmospheric model calculations 
(Grobman and Ingebo, 1974) . 

An upper limit to the emission index may be evaluated 
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based on stoichiometric considerations assuming perfect 
combustion efficiency. For example, for a kerosene mixture 
(Jet A) , a water vapor emission index EI(H20) of about 
1400 g/kg fuel is predicted. Higher order distillates such 
as naptha produce slightly higher yields of H 2 O*. Measure- 
ments indicate slightly lower values of emission index 
(CIAP, 1975) and 1250 g/kg fuel is the recommended value 
for cruise conditions and all engine types. 

The emission of sulfur, primarily as SC^/ is deter- 
mined by the fuel stock sulfur content which is expressed 
as percent by weight, % (w/w) . This abundance is controlla- 
ble by implementing available fuel desulfurization techni- 
ques. Current average aviation fuel sulfur contents are 
0.05 percent (w/w) with a maximum allowable value of 0.30 
percent (w/w) (CIAP, 1975) . These values suggest typical 
and worse case EI(S) of 0.5 and 3.0 g/kg fuel, respectively. 

Elevated combustion zone temperatures and the presence 
of radicals such as O, N, and OH facilitate the formation 
of oxides of nitrogen which would not be energetically 
favored under ambient conditions. Since NO x emissions are 
highly dependent on operating conditions , values at cruise 
conditions are assumed. It is necessary to distinguish 
between near term (present day) and far term design con- 
cepts in assigning NO x emission indices (Poppoff et al. , 
1978) since future technological advances are planned to 
reduce N0 X emissions by a factor of 6 (Reck, 1978) . 


* 

It is notable that H 2 O emission indices for coal-derived 
fuels suggested for implementation in certain advanced SST 
(ASST) design concepts (Witcofski, 1978) are significantly 
higher. For example, synthetic liquid natural gas (SYN-LNG) 
has a El (H 2 O) of about 2300 g/kg fuel while that for liquid 
hydrogen (LH 2 ) is near 9000 g/kg fuel. On an energy equiv- 
alent basis, these emissions are 60 and 160 percent higher, 
respectively, than Jet A (Witcofski, 1978). 
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Typical near term emissions [EI(NO x ) = 18 g/kg fuel] are 
based on the Olympus 593 engine as described in Lohmann and 
Riecke (1977) ; this value is consistent with that given by 
Grobman and Ingebo (1974) . The corresponding far term 
technology estimate is 3 g/kg fuel. Values of El (N0 X ) are 
summarized in Table F-l. 

Assessment of potential perturbations to atmospheric 
chemistry and global climate due to aircraft exhaust 
emissions is dependent on the total fleet-derived source 
strength, R^ (MT yr ^) for each pollutant species, which is 
defined as 

R ± - (El) :L M$T f /10' L2 (F-l) 

where M is the number of aircraft in the proposed nominal 
fleet, $ (kg fuel hr“l) the engine fuel flow rate at cruise 
altitude and T^ (hr/yr) the yearly cruise time in hours for 
each aircraft. Typical source strengths based on the given 
emission indices and a nominal fleet of 500 aircraft are 
evaluated and presented in Table F-l. In these calcula- 
tions, we assume an average engine fuel flow rate of 
3.78x10^ kg fuel hr”-*-, which is based on a "Type A" air- 
craft as defined in Poppoff et al. (1978) and described in 
Baber and Swanson (1976) . Cruise time is taken to be 7 hr 
per day per aircraft. 

Equivalent aircraft production rates P a are evaluated 
assuming injection of N0 X as pure NO^ and distributing it 
uniformly over the globe in a 1 km thick spherical shell 
(to be centered at the cruise altitudes of 15.2 and 


+ • 

The assumption of pure NO production is only approximately 

true. The volume fraction of NO, however, is equal to 

90-95 percent (CIAP, 1975) of the total NO x emissions. 
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18.3 km). Based on these assumptions, the given aircraft 
characteristics are converted to equivalent production 
rates using: 

P a = (7 . 82 x10“ 2 )M(EI) (F-2 ) 

Calculations based on equations (F-l) and F-2) for N0 X are 
summarized in Table F-2. Ozone perturbation sensitivity 
analysis will be performed on these values. 
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